假设有必要计算打包浮点数据的倒数或倒数平方根。两者都可以轻松完成：
__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.5 * 2 ^ -12。鉴于单精度浮点数的机器epsilon是2 ^ -24，我们可以说这个近似值近似 半 精确。
似乎可以添加Newton-Raphson迭代来产生结果 单 精度（可能不如IEEE标准所要求的那样精确），请参阅 GCC， ICC，讨论 LLVM。从理论上讲，同样的方法可以用于双精度值，产生 半 要么 单 要么 双 精确。
我有兴趣为float和double数据类型以及所有（half，single，double）精度实现此方法的实现。处理特殊情况（除以零，sqrt（-1），inf / nan等）不是必需的。此外，我不清楚这些例程中哪一个比普通的IEEE编译解决方案更快，哪个更慢。
以下是对答案的一些小限制，请：
- 在代码示例中使用内在函数。程序集依赖于编译器，因此不太有用。
- 对函数使用类似的命名约定。
- 实现例程，将包含密集打包的float / double值的单个SSE / AVX寄存器作为输入。如果有相当大的性能提升，你也可以发布几个寄存器作为输入的例程（两个reg可能是可行的）。
- 如果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迭代。
- RECIP 上 浮动需要 1,4 周期与 7 周期。
- rsqrt 上 浮动需要 1,6 周期与 14 周期。
- RECIP 上 双需要 3,6,9 周期与 14 周期。
- rsqrt 上 双需要 3,8,13 周期与 28 周期。
警告：我必须创造性地将原始结果舍入...
在实践中有很多算法的例子。例如：
Newton Raphson与SSE2 - 有人可以解释我这3行 有一个答案解释使用的迭代 英特尔的一个例子。
对于井上分析，让我们说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.
如果它不是依赖链的一部分，那么在这里有很多其他计算重叠的空间。但是，如果代码的每次迭代的数据都依赖于前面的数据，那么11周期的延迟可能会更好。 sqrtps。或者即使每次循环迭代足够长，乱序执行也无法通过重叠独立迭代来隐藏它们。
要得到 sqrt(x) 代替 1/sqrt(x)， 乘以 x （并修复如果 x 可以是零，例如通过掩蔽（_mm_andn_ps）使用CMPPS对0.0的结果。最佳方式是更换 half_nr 同 half_xnr = _mm_mul_ps( half, xnr );。这并没有延长部门链，因为 half_xnr 可以在第11周期开始，但直到结束（第19周期）才需要。与可用的FMA相同：总延迟没有增加。
如果11位精度足够（没有牛顿迭代）， 英特尔优化手册（第11.12.3节） 建议使用rcpps（rsqrt（x）），它比原来的x更差，至少使用AVX。它可能会保存带有128位SSE的movdqa指令，但256b rcpps比256b mul或fma慢。 （并且它允许您在最后一步使用FMA而不是mul免费添加sqrt结果。）
这个循环的SSE版本，不考虑任何移动指令，是6微秒。这意味着它应该具有Haswell的吞吐量 每3个循环一个 （假设两个执行端口可以处理FP mul，而rsqrt在FP add / sub的相对端口上）。在SnB / IvB（可能还有Nehalem）上，它的吞吐量应该是 每5个循环一个，因为mulps和rsqrtps竞争端口0.subps在port1上，这不是瓶颈。
对于Haswell，我们可以使用FMA将减法与mul结合起来。但是，分频器/ sqrt单元不是256b宽，因此与其他所有内容不同，ymm regs上的divps / sqrtps / rsqrtps / rcpps需要额外的uops和额外的周期。
// 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-cyle-slow 256b rsqrt来抵消，净延迟减少1个周期（相当于两倍宽）。 SnB / IvB AVX将是128b的24c延迟，256b的26c延迟。
256b FMA版本总共使用7个uop。 （VRSQRTPS是2 uop，2是p0，一个是p1 / 5.）256b mulps和fma都是单uop指令，两者都可以在端口0或端口1上运行。（p0仅在Haswell之前）。吞吐量应该是： 每3c一个，如果OOO引擎将uop调度到最佳执行端口。 （即来自rsqrt的shuffle uop总是进入p5，永远不会进入p1，它将占用mul / fma带宽。）至于与其他计算重叠（不仅仅是自身的独立执行），它再次非常轻量级。只有7个具有23个周期的dep链的uops为其他东西留下了很大的空间，而那些uops位于重新排序的缓冲区中。
如果这是一个巨大的dep链中的一步，没有其他任何事情发生（甚至不是一个独立的下一次迭代），那么256b vsqrtps 是19周期延迟，每14个周期一个吞吐量。 （Haswell的）。如果你仍然真的需要倒数，那么256b vdivps 也有18-21c的延迟，每14c吞吐量一个。因此对于普通sqrt，指令具有较低的延迟。对于recip sqrt，近似的迭代是较少的延迟。 （如果它与自身重叠，那么FAR会增加吞吐量。如果与其他没有分隔单元的东西重叠， sqrtps 不是问题。）
sqrtps 可以是吞吐量赢 rsqrt +牛顿迭代，如果它是一个循环体的一部分，其中有足够的其他非分割和非sqrt工作，则除法单元不饱和。
如果您需要，尤其如此 sqrt(x)不是 1/sqrt(x)。例如在Haswell上使用AVX2，在一个适合L3缓存的浮点数组上的复制+ arcsinh循环，实现为 fastlog(v + sqrt(v*v + 1))，使用真实VSQRTPS或使用VRSQRTPS + Newton-Raphson迭代以大约相同的吞吐量运行。 （即使有一个 log（）的非常快速的近似值，所以总循环体大约是9 FMA / add / mul / convert操作，2 boolean，再加上VSQRTPS ymm。使用刚才有一个加速 fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1)因此它不会对内存带宽产生瓶颈，但它可能会导致延迟出现瓶颈（因此无序执行无法利用独立迭代的所有并行性）。
其他精度
对于 半精度，没有指示在半浮点数上进行计算。您打算在使用转换说明加载/存储时动态转换。
对于 双精度，没有 rsqrtpd。据推测，我们的想法是，如果你不需要完全精确，你应该首先使用浮动。所以你可以转换为float和back，然后执行完全相同的算法，但是 pd 代替 ps 内部函数。或者，您可以将数据保持浮动一段时间。例如将两个ymm的双重寄存器转换成一个单独的ymm寄存器。
不幸的是，没有一条指令需要两个ymm的双精度寄存器并输出一个单独的ymm。那你必须去ymm-> xmm两次 _mm256_insertf128_ps 一个xmm到另一个的128。但是，您可以将256b ymm向量提供给相同的算法。
如果你要在之后转换回加倍，那么分别对两个单独的128b寄存器进行rsqrt + Newton-Raphson迭代是有意义的。插入/提取的额外2个uop，以及256b rsqrt的额外2个uop，开始加起来，更不用说3个周期的延迟了 vinsertf128 / vextractf128。计算将重叠，两个结果将分开几个周期。 （或者相隔1个周期，如果您有一个特殊版本的代码，可以同时对2个输入进行交错操作）。
请记住，单精度的指数范围小于double，因此转换可以溢出到+ Inf或下溢到0.0。如果你不确定，一定只是使用正常 _mm_sqrt_pd。
有了AVX512F，就有了 _mm512_rsqrt14_pd( __m512d a)。同 AVX512ER（KNL但不是SKX或Cannonlake）， 有 _mm512_rsqrt28_pd。 _ps 版本当然也存在。在更多情况下，14位尾数精度可能足以在没有牛顿迭代的情况下使用。
牛顿迭代仍然不会给你一个 正确舍入 单精度浮点数就像常规sqrt一样。
在实践中有很多算法的例子。例如：
Newton Raphson与SSE2 - 有人可以解释我这3行 有一个答案解释使用的迭代 英特尔的一个例子。
对于井上分析，让我们说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.
如果它不是依赖链的一部分，那么在这里有很多其他计算重叠的空间。但是，如果代码的每次迭代的数据都依赖于前面的数据，那么11周期的延迟可能会更好。 sqrtps。或者即使每次循环迭代足够长，乱序执行也无法通过重叠独立迭代来隐藏它们。
要得到 sqrt(x) 代替 1/sqrt(x)， 乘以 x （并修复如果 x 可以是零，例如通过掩蔽（_mm_andn_ps）使用CMPPS对0.0的结果。最佳方式是更换 half_nr 同 half_xnr = _mm_mul_ps( half, xnr );。这并没有延长部门链，因为 half_xnr 可以在第11周期开始，但直到结束（第19周期）才需要。与可用的FMA相同：总延迟没有增加。
如果11位精度足够（没有牛顿迭代）， 英特尔优化手册（第11.12.3节） 建议使用rcpps（rsqrt（x）），它比原来的x更差，至少使用AVX。它可能会保存带有128位SSE的movdqa指令，但256b rcpps比256b mul或fma慢。 （并且它允许您在最后一步使用FMA而不是mul免费添加sqrt结果。）
这个循环的SSE版本，不考虑任何移动指令，是6微秒。这意味着它应该具有Haswell的吞吐量 每3个循环一个 （假设两个执行端口可以处理FP mul，而rsqrt在FP add / sub的相对端口上）。在SnB / IvB（可能还有Nehalem）上，它的吞吐量应该是 每5个循环一个，因为mulps和rsqrtps竞争端口0.subps在port1上，这不是瓶颈。
对于Haswell，我们可以使用FMA将减法与mul结合起来。但是，分频器/ sqrt单元不是256b宽，因此与其他所有内容不同，ymm regs上的divps / sqrtps / rsqrtps / rcpps需要额外的uops和额外的周期。
// 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-cyle-slow 256b rsqrt来抵消，净延迟减少1个周期（相当于两倍宽）。 SnB / IvB AVX将是128b的24c延迟，256b的26c延迟。
256b FMA版本总共使用7个uop。 （VRSQRTPS是2 uop，2是p0，一个是p1 / 5.）256b mulps和fma都是单uop指令，两者都可以在端口0或端口1上运行。（p0仅在Haswell之前）。吞吐量应该是： 每3c一个，如果OOO引擎将uop调度到最佳执行端口。 （即来自rsqrt的shuffle uop总是进入p5，永远不会进入p1，它将占用mul / fma带宽。）至于与其他计算重叠（不仅仅是自身的独立执行），它再次非常轻量级。只有7个具有23个周期的dep链的uops为其他东西留下了很大的空间，而那些uops位于重新排序的缓冲区中。
如果这是一个巨大的dep链中的一步，没有其他任何事情发生（甚至不是一个独立的下一次迭代），那么256b vsqrtps 是19周期延迟，每14个周期一个吞吐量。 （Haswell的）。如果你仍然真的需要倒数，那么256b vdivps 也有18-21c的延迟，每14c吞吐量一个。因此对于普通sqrt，指令具有较低的延迟。对于recip sqrt，近似的迭代是较少的延迟。 （如果它与自身重叠，那么FAR会增加吞吐量。如果与其他没有分隔单元的东西重叠， sqrtps 不是问题。）
sqrtps 可以是吞吐量赢 rsqrt +牛顿迭代，如果它是一个循环体的一部分，其中有足够的其他非分割和非sqrt工作，则除法单元不饱和。
如果您需要，尤其如此 sqrt(x)不是 1/sqrt(x)。例如在Haswell上使用AVX2，在一个适合L3缓存的浮点数组上的复制+ arcsinh循环，实现为 fastlog(v + sqrt(v*v + 1))，使用真实VSQRTPS或使用VRSQRTPS + Newton-Raphson迭代以大约相同的吞吐量运行。 （即使有一个 log（）的非常快速的近似值，所以总循环体大约是9 FMA / add / mul / convert操作，2 boolean，再加上VSQRTPS ymm。使用刚才有一个加速 fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1)因此它不会对内存带宽产生瓶颈，但它可能会导致延迟出现瓶颈（因此无序执行无法利用独立迭代的所有并行性）。
其他精度
对于 半精度，没有指示在半浮点数上进行计算。您打算在使用转换说明加载/存储时动态转换。
对于 双精度，没有 rsqrtpd。据推测，我们的想法是，如果你不需要完全精确，你应该首先使用浮动。所以你可以转换为float和back，然后执行完全相同的算法，但是 pd 代替 ps 内部函数。或者，您可以将数据保持浮动一段时间。例如将两个ymm的双重寄存器转换成一个单独的ymm寄存器。
不幸的是，没有一条指令需要两个ymm的双精度寄存器并输出一个单独的ymm。那你必须去ymm-> xmm两次 _mm256_insertf128_ps 一个xmm到另一个的128。但是，您可以将256b ymm向量提供给相同的算法。
如果你要在之后转换回加倍，那么分别对两个单独的128b寄存器进行rsqrt + Newton-Raphson迭代是有意义的。插入/提取的额外2个uop，以及256b rsqrt的额外2个uop，开始加起来，更不用说3个周期的延迟了 vinsertf128 / vextractf128。计算将重叠，两个结果将分开几个周期。 （或者相隔1个周期，如果您有一个特殊版本的代码，可以同时对2个输入进行交错操作）。
请记住，单精度的指数范围小于double，因此转换可以溢出到+ Inf或下溢到0.0。如果你不确定，一定只是使用正常 _mm_sqrt_pd。
有了AVX512F，就有了 _mm512_rsqrt14_pd( __m512d a)。同 AVX512ER（KNL但不是SKX或Cannonlake）， 有 _mm512_rsqrt28_pd。 _ps 版本当然也存在。在更多情况下，14位尾数精度可能足以在没有牛顿迭代的情况下使用。
牛顿迭代仍然不会给你一个 正确舍入 单精度浮点数就像常规sqrt一样。