根据精度，使用SSE/AVX快速矢量化的rsqrt和倒数。

Question

根据精度，使用SSE/AVX快速矢量化的rsqrt和倒数。

16

假设需要计算打包的浮点数据的倒数或倒数平方根。两者都可以通过以下方式轻松完成：

__m128 recip_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), x); }
__m128 rsqrt_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(x)); }

这个代码可以正常工作但速度较慢：根据指南，在Sandy Bridge上需要14和28个周期（吞吐量）。相应的AVX版本在Haswell上需要几乎相同的时间。

另一方面，可以使用以下版本代替：

__m128 recip_float4_half(__m128 x) { return _mm_rcp_ps(x); }
__m128 rsqrt_float4_half(__m128 x) { return _mm_rsqrt_ps(x); }

他们只需要1到2个周期的时间（吞吐量），从而大大提高性能。然而，它们非常近似：它们产生的结果相对误差小于1.5 * 2^-12。考虑到单精度浮点数的机器epsilon为2^−24，我们可以说这种近似具有约半精度。

看起来可以添加牛顿-拉弗森迭代以产生带有单精度的结果（可能不如IEEE标准要求的那么精确），请参见GCC，ICC，以及LLVM上的讨论。理论上，相同的方法可以用于双精度值，产生半精度或单精度或双精度。

我对此方法的实现很感兴趣，包括浮点和双精度数据类型以及所有（半精度、单精度、双精度）精度。处理特殊情况（除零、sqrt（-1）、inf / nan等）并不是必要的。此外，我不清楚这些例程中哪些比微不足道的IEEE兼容解决方案更快，哪些会更慢。

以下是答案的一些小限制，请注意：

在代码示例中使用内部函数。汇编语言依赖于编译器，因此不太有用。
使用类似的函数命名约定。
实现以单个SSE / AVX寄存器为输入的紧密打包的float / double值的例程。如果有相当大的性能提升，您还可以发布采用几个寄存器作为输入的例程（两个寄存器可能可行）。
如果它们绝对相等，不要同时发布SSE / AVX版本，直到更改_mm为_mm256或反之亦然。

欢迎任何性能估计、测量和讨论。

总结：

以下是单精度浮点数的版本，其中包含一个NR迭代：

__m128 recip_float4_single(__m128 x) {
  __m128 res = _mm_rcp_ps(x);
  __m128 muls = _mm_mul_ps(x, _mm_mul_ps(res, res));
  return res =  _mm_sub_ps(_mm_add_ps(res, res), muls);
}
__m128 rsqrt_float4_single(__m128 x) {
  __m128 three = _mm_set1_ps(3.0f), half = _mm_set1_ps(0.5f);
  __m128 res = _mm_rsqrt_ps(x); 
  __m128 muls = _mm_mul_ps(_mm_mul_ps(x, res), res); 
  return res = _mm_mul_ps(_mm_mul_ps(half, res), _mm_sub_ps(three, muls)); 
}

Peter Cordes提供的答案解释了如何创建其他版本，并包含了彻底的理论性能分析。

您可以在这里找到所有实现的解决方案和基准测试结果：recip_rsqrt_benchmark。

在Ivy Bridge上获得的吞吐量结果如下所示。仅对单寄存器SSE实现进行了基准测试。时间以每次调用的周期数给出。第一个数字是半精度（无NR），第二个数字是单精度（1个NR迭代），第三个数字是2个NR迭代。

float 上的 recip 需要 1, 4 个周期，而不是 7 个周期。
float 上的 rsqrt 需要 1, 6 个周期，而不是 14 个周期。
double 上的 recip 需要 3, 6, 9 个周期，而不是 14 个周期。
double 上的 rsqrt 需要 3, 8, 13 个周期，而不是 28 个周期。

警告：我必须创造性地四舍五入原始结果...

- stgatilov

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Peter Cordes · Accepted Answer

有许多实际应用算法的例子。例如：

牛顿-拉弗森算法与SSE2 - 有人能解释一下这三行吗？提供了一个解释Intel示例之一所使用的迭代。

对于性能分析，以Haswell为例（它可以在两个执行端口上进行FP乘法，而不像以前的设计），我将在此重现代码（每行一个操作符）。请参见http://agner.org/optimize/查看指令吞吐量和延迟表，并了解更多背景信息的文档。

// execute (aka dispatch) on cycle 1, results ready on cycle 6
nr = _mm_rsqrt_ps( x );

// both of these execute on cycle 6, results ready on cycle 11
xnr = _mm_mul_ps( x, nr );         // dep on nr
half_nr = _mm_mul_ps( half, nr );  // dep on nr

// can execute on cycle 11, result ready on cycle 16
muls = _mm_mul_ps( xnr , nr );     // dep on xnr

// can execute on cycle 16, result ready on cycle 19
three_minus_muls = _mm_sub_ps( three, muls );  // dep on muls

// can execute on cycle 19, result ready on cycle 24
result = _mm_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls
// result is an approximation of 1/sqrt(x), with ~22 to 23 bits of precision in the mantissa.

如果其他计算与该处的依赖关系无关，则有很大的重叠空间。但是，如果您代码中每次迭代的数据依赖于先前的数据，则使用sqrtps的11个周期延迟可能更好。或者即使每个循环迭代都足够长，使得乱序执行无法通过重叠独立迭代来隐藏所有操作，也可以选择使用sqrtps。

为了得到sqrt(x)而不是1/sqrt(x)，请乘以x（如果x可能为零，则需要修复，例如通过使用CMPPS和0.0进行屏蔽(_mm_andn_ps)的结果）。最优方式是用half_xnr = _mm_mul_ps( half, xnr );替换half_nr。这不会增加依赖链的长度，因为half_xnr可以从第11个周期开始，但直到结束（第19个周期）才需要。如果FMA可用，则情况相同：总延迟不增加。

如果11位精度足够（无牛顿迭代），Intel的优化手册（第11.12.3节）建议使用rcpps(rsqrt(x))，但至少在AVX下，这比乘以原始x更差。它可能会在128位SSE中节省一个movdqa指令，但256b rcpps比256b mul或fma慢。（并且它让你可以将sqrt结果免费添加到FMA中而不是mul作为最后一步）。

该循环的 SSE 版本，不考虑任何移动指令，需要 6 个微操作。这意味着在 Haswell 上它的吞吐量应该是 每 3 个周期完成一次（因为两个执行端口可以处理 FP mul，rsqrt 在 FP add/sub 的相反端口）。在 SnB/IvB（和可能是 Nehalem）上，由于 mulps 和 rsqrtps 竞争使用端口0，其吞吐量应该为每 5 个周期完成一次。subps 在端口1上，不是瓶颈。

对于 Haswell，我们可以使用 FMA 将减法与乘法结合起来。然而，除法器/平方根单元不是 256b 宽，因此与其他所有内容不同，在 ymm 寄存器上进行的 divps / sqrtps / rsqrtps / rcpps 需要额外的微操作和额外的周期。

// vrsqrtps ymm has higher latency
// execute on cycle 1, results ready on cycle 8
nr = _mm256_rsqrt_ps( x );

// both of can execute on cycle 8, results ready on cycle 13
xnr = _mm256_mul_ps( x, nr );         // dep on nr
half_nr = _mm256_mul_ps( half, nr );  // dep on nr

// can execute on cycle 13, result ready on cycle 18
three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three );  // -(xnr*nr) + 3

// can execute on cycle 18, result ready on cycle 23
result = _mm256_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls

我们使用FMA可以节省3个周期。虽然使用了2个周期较慢的256b rsqrt，但净收益是少1个周期的延迟（对于两倍宽度来说相当不错）。SnB/IvB AVX在128b时延迟为24c，在256b时延迟为26c。

256b FMA版本总共使用7个微操作。（VRSQRTPS是3个微操作，2个用于p0，1个用于p1/5。）256b mulps和fma都是单微操作指令，都可以在端口0或端口1上运行。（仅在Haswell之前的p0上。）因此，吞吐量应该是：每3个周期一个，如果OOO引擎将微操作分派到最佳执行端口。（即rsqrt的shuffle微操作始终进入p5，而不是占用mul/fma带宽的p1。）至于与其他计算重叠（不仅是独立执行自身），再次说明它非常轻量级。只有7个微操作和23个周期的依赖链，在重新排序缓冲区中等待这些微操作时留下了很多空间供其他事情发生。

如果这是巨大依赖链中的一步，除了独立的下一次迭代之外没有其他事情发生，那么256b的vsqrtps的延迟为19个周期，吞吐量为每14个周期一个（Haswell）。如果您仍然需要倒数，则256b的vdivps的延迟也为18-21c，吞吐量为每14c一个。因此，对于普通的平方根，指令具有较低的延迟。对于倒数平方根，近似的迭代具有更少的延迟。（如果与自身重叠，则吞吐量要高得多。如果与不使用除法单元的其他内容重叠，则sqrtps不是问题。）

如果sqrtps是循环体的一部分，并且有足够的其他非除法和非平方根工作正在进行，使除法单元未饱和，则sqrtps可以成为吞吐量优势，而不是rsqrt+牛顿迭代。

这尤其适用于需要sqrt(x)而不是1/sqrt(x)的情况。例如，在Haswell with AVX2上，对于一个适合L3缓存的浮点数数组的复制+arcsinh循环，实现为fastlog(v + sqrt(v*v + 1))，使用真实的VSQRTPS或VRSQRTPS + Newton-Raphson迭代，吞吐量大致相同。(即使使用非常快的对数近似,因此总循环体大约有9个FMA/add/mul/convert操作和2个布尔操作，加上VSQRTPS ymm。只使用fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1)会有加速，所以它不会受到内存带宽的瓶颈限制，但它可能会受到延迟的瓶颈限制(因此乱序执行无法利用独立迭代的所有并行性)。

其他精度

对于半精度，没有指令可以对半浮点数进行计算。您需要在加载/存储时使用转换指令进行即时转换。

对于双精度，没有rsqrtpd。大概是因为如果您不需要完全精度，那么应该首先使用float。所以您可以将其转换为float，然后再次使用相同的算法，但使用pd而不是ps内部函数。或者，您可以将数据保留为float一段时间。例如，将两个双倍长度的ymm寄存器转换为一个单倍长度的ymm寄存器。

不幸的是，没有单个指令可以将两个双倍长度的ymm寄存器输出为一个单倍长度的ymm。您需要两次将ymm->xmm，然后将一个xmm插入到另一个xmm的高128位中。但然后您可以将该256b ymm向量提供给相同的算法。

如果您在之后会转换为double类型，那么将rsqrt和牛顿-拉夫逊迭代分别在两个128b的单精度寄存器上进行可能更加合理。插入/提取的额外2个uops以及256b rsqrt的额外2个uops开始累加，更不用说vinsertf128 / vextractf128的3个周期延迟。计算将重叠，并且两个结果将准备好一些循环之后（如果您有一个特殊版本的代码来交错同时对两个输入执行操作，则为1个循环之后）。

请记住，单精度具有比double小的指数范围，因此转换可能会溢出到+Inf或下溢到0.0。如果您不确定，请使用正常的_mm_sqrt_pd。

使用AVX512F，有_mm512_rsqrt14_pd(__m512d a)。使用AVX512ER（仅限KNL而不是SKX或Cannonlake），有_mm512_rsqrt28_pd。当然也存在_ps版本。在更多情况下，14位尾数精度可能足以在没有牛顿迭代的情况下使用。

牛顿迭代仍无法像常规sqrt一样为您提供正确舍入的单精度浮点数。