快速向量化rsqrt和SSE / AVX取决于精度的倒数
假设有必要计算打包浮点数据的倒数或倒数平方根。 两者可以通过以下方式轻松完成
__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。 鉴于单精度浮点数的机器ε为2 ^ -24，我们可以说这个近似具有大约一半的精度。
似乎可以添加Newton-Raphson迭代来产生单精度结果（可能不像IEEE标准所要求的那样精确），参见GCC和ICC在LLVM上的讨论。 理论上，对于双精度值可以使用相同的方法，产生一半或单精度或双精度。
我对float和double数据类型以及所有（半，单，双）精度的这种方法的工作实现感兴趣。 处理特殊情况（除零，sqrt（-1），inf / nan等）是不必要的。 另外，对于我来说，这些例程中的哪一个比简单的IEEE编译器解决方案更快，哪个更慢，这一点并不清楚。
以下是对答案的一些小的限制，请：
任何绩效评估，测量和讨论都是受欢迎的。
概要
下面是一个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迭代。
警告：我不得不创造性地整理原始结果...
实践中有很多算法的例子。 例如：
Newton Raphson与SSE2--有人能解释我这三行有一个答案，解释了英特尔的一个例子所使用的迭代。
为了进行性能分析，我们假设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个周期延迟可能会更好。 或者即使每个循环迭代足够长以至于乱序执行无法通过重叠独立迭代来隐藏全部循环。
为了得到sqrt(x)而不是1/sqrt(x) ，乘以x （并且如果x可以为0，例如通过掩蔽（ _mm_andn_ps ）并且_mm_andn_ps的结果为0.0）。 最佳方法是用half_xnr = _mm_mul_ps( half, xnr );替换half_nr half_xnr = _mm_mul_ps( half, xnr ); 。 这不会延长dep链，因为half_xnr可以从第11个循环开始，但直到结束（循环19）才需要。 与FMA可用相同：总延迟没有增加。
如果11位精度足够高（无牛顿迭代），英特尔优化手册（第11.12.3节）建议使用rcpps（rsqrt（x）），这比至少与AVX相乘的结果要差得多。 它可以保存128位SSE的movdqa指令，但256b rcpps比256b mul或fma慢。 （并且它可以让您将sqrt结果添加到FMA的免费项目中，而不是最后一步的mul项目中）。
这个循环的SSE版本不包含任何移动指令，是6个uops。 这意味着它应该在Haswell上的吞吐量为每3个周期一个 （因为两个执行端口可以处理FP mul，并且rsqrt位于与FP add / sub相反的端口上）。 在SnB / IvB（可能是Nehalem）上， 每5个周期应该有一个吞吐量，因为mulps和rsqrtps竞争端口0. subps在端口1上，这不是瓶颈。
对于Haswell，我们可以使用FMA将减法与mul结合。 然而，分频器/ sqrt单元不是256b宽，因此不同于其他的，ymm regs上的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-cyle-256b rsqrt来弥补，因为净增益为1周期更少的延迟（相当好，是两倍宽）。 SnB / IvB AVX将是128b的24c延迟，256b的26c延迟。
256b FMA版本总共使用了7个uops。 （ VRSQRTPS为3 VRSQRTPS ，2为p0，另一个为p1 / 5）256b mulps和fma都是单端指令，都可以在端口0或端口1上运行（仅在前Haswell上为p0）。 因此吞吐量应该是：如果OOO引擎将uops分派到最佳执行端口，则每3c一个吞吐量。 （也就是说，来自rsqrt的shuffle uop总是进入p5，而不是p1，它将占用多个/ fma带宽。）至于与其他计算重叠（不仅仅是独立执行自身），它又是非常轻量级的。 只有7架uops带有23个循环的dep链，在这些uoop坐在重新排序缓冲区中时，会留下很多空间让其他东西发生。
如果这是一个巨型dep链中的一步，没有其他任何事情发生（甚至不是独立的下一次迭代），那么256b vsqrtps是19周期延迟，吞吐量为每14个周期一个。 （Haswell的）。 如果你真的需要互惠，那么256b vdivps也有18-21c的延迟，每14c吞吐量有一个。 所以对于普通的sqrt，指令的延迟较低。 对于收信人sqrt，迭代的近似是较少的延迟。 （如果它与自身重叠，那么FAR会有更多的吞吐量，如果与其他不是分频单元的东西重叠， sqrtps不是问题。）
sqrtps可以是吞吐量胜利vs. rsqrt +牛顿迭代，如果它是循环主体的一部分，并且有足够的其他非分割和非sqrt工作，分割单元不会饱和。
如果您需要sqrt(x) ，而不是1/sqrt(x) ，则尤其如此 。 例如在具有AVX2的Haswell上，通过实施为fastlog(v + sqrt(v*v + 1))适合L3缓存的浮点数组上的副本+ arcsinh循环，以与真实VSQRTPS或VRSQRTPS +牛顿 - 拉夫逊迭代。 （即使对log（）有非常快速的近似，所以总循环体大约为9 FMA / add / mul / convert操作，2布尔值加上VSQRTPS ymm。使用fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1) ，所以它不是内存带宽的瓶颈，但它可能是延迟瓶颈（因此无序执行不能利用独立迭代的所有并行）。
其他精度
对于半精度 ，没有指令在半浮点上进行计算。 您打算在加载/存储时使用转换说明即时转换。
对于双精度 ，没有rsqrtpd 。 据推测，这个想法是，如果你不需要完全精确的话，你应该首先使用float。 所以你可以将其转换为float和back，然后执行完全相同的算法，但使用pd而不是ps内在函数。 或者，您可以将数据保持浮动一段时间。 例如将两个ymm的双打寄存器转换成一个ymm的单打寄存器。
不幸的是，没有一条指令需要两个ymm的双打寄存器并输出一个ymm的单数。 你必须两次走ymm-> xmm，然后_mm256_insertf128_ps一个xmm到另一个的最高128位。 但是，您可以将该256b ymm矢量提供给相同的算法。
如果你打算在之后立即转换为double，那么分别对两个128b寄存器单独进行rsqrt + Newton-Raphson迭代可能是有意义的。 额外的2个uops插入/提取，额外的2个uops用于256b rsqrt，开始累计，更不用说vinsertf128 / vextractf128的3个周期延迟。 计算会重叠，两个结果将准备好几个周期。 （或者，如果您有一个特殊版本的代码用于一次输入2个输入的交叉操作）
请记住，单精度指数的范围小于双精度，因此转换可以溢出到+ Inf或下溢到0.0。 如果您不确定，那么只需使用正常的_mm_sqrt_pd 。
使用AVX512F时，有_mm512_rsqrt14_pd( __m512d a) 。 使用AVX512ER（KNL而不是SKX或Cannonlake），有_mm512_rsqrt28_pd 。 _ps版本当然也存在。 在更多情况下，14位的尾数精度可能足以在没有牛顿迭代的情况下使用。
牛顿迭代仍然不会给你一个正常的四舍五入单精度浮点的方式常规sqrt会。
链接地址: http://www.djcxy.com/p/85643.html上一篇: Fast vectorized rsqrt and reciprocal with SSE/AVX depending on precision
下一篇: Floating point division vs floating point multiplication