在_mm256_rsqrt_ps（）中处理零
鉴于_mm256_sqrt_ps()相对较慢，并且我生成的值会立即被_mm256_floor_ps()截断，看看它似乎在做：
_mm256_mul_ps(_mm256_rsqrt_ps(eightFloats), eightFloats);
是获得额外性能和避免管道停滞的方法。
不幸的是，零值，我当然得到崩溃计算1/sqrt(0) 。 围绕这个的最佳方法是什么？ 我试过这个（有效并且更快），但是有更好的方法，还是我会在某些条件下遇到问题？
_mm256_mul_ps(_mm256_rsqrt_ps(_mm256_max_ps(eightFloats, _mm256_set1_ps(0.1))), eightFloats);
我的代码用于垂直应用程序，因此我可以假设它将在Haswell CPU（i7-4810MQ）上运行，因此可以使用FMA / AVX2。 原始代码大约是：
float vals[MAX]; int sum = 0; for (int i = 0; i < MAX; i++) { int thisSqrt = (int) floor(sqrt(vals[i])); sum += min(thisSqrt, 0x3F); }
vals所有值都应该是整数值。 （为什么一切都不仅仅是int是一个不同的问题…）
tl; dr ：查看编译和应该工作的代码的结尾。
要解决0.0问题，您还可以使用0.0特殊情况输入，其中源与FP的FP比较。 使用比较结果作为掩码，将sqrt(x) = x * rsqrt(x) 0 * +Infinity任何NaN归零。 Clang在自动向量化时执行此操作。 （但它使用具有归零向量的混合，而不是使用带有andnps的比较掩码直接归零或保留元素。）
如njuffa所建议的那样 ，也可以使用sqrt(x) ~= recip(rsqrt(x)) 。 rsqrt(0) = +Inf 。 recip(+Inf) = 0 。 但是，使用两个近似值会使相对误差复杂化，这是一个问题。
你缺少的东西：
当输入是完美的正方形时，截断为整数（而不是舍入）需要精确的sqrt结果。 如果25 * rsqrt（25）的结果是4.999999或其他东西（而不是5.00001），你将添加4而不是5 。
即使使用Newton-Raphson迭代， rsqrtps 也不像 rsqrtps 那样完全准确 ，因此它仍然可以提供5.0 - 1ulp 。 （1ulp =最后一位的一个单位=尾数的最低位）。
也：
- Newton Raphson公式解释道
- Newton Raphson SSE实施性能（延迟/吞吐量） 。 请注意，我们更关心吞吐量而不是延迟，因为我们在一个没有做太多其他事情的循环中使用它。 sqrt不是循环承载的dep链的一部分，因此不同的迭代可以同时在其中运行sqrt计算。
在进行(x+offset)*approx_rsqrt(x+offset)然后截断为整数之前，可以通过添加一个小常量来杀死2只一石。 大到足以克服1.5 * 2 -12的最大相对误差，但小到足以不使sqrt_approx(63*63-1+offset)高达63 （最敏感的情况）。
63*1.5*2^(-12) == 0.023071... approx_sqrt(63*63-1) == 62.99206... +/- 0.023068..
实际上，即使没有添加任何内容，我们也会在没有牛顿迭代的情况下搞砸 。 approx_sqrt(63*63-1)可以approx_sqrt(63*63-1)出现在63.0以上。 n=36是sqrt(n*n-1) + error中的相对误差小于sqrt(n*n)的最大值。 GNU Calc：
define f(n) { local x=sqrt(n*n-1); local e=x*1.5*2^(-12); print x; print e, x+e; } ; f(36) 35.98610843089316319413 ~0.01317850650545403926 ~35.99928693739861723339 ; f(37) 36.9864840178138587015 ~0.01354485498699237990 ~37.00002887280085108140
您的源数据是否具有任何属性，这意味着您不必担心它只是在一个大的完美正方形下面？ 例如它总是完美的方块吗？
您可以检查所有可能的输入值，因为重要的域非常小（整数FP值从0..63 * 63），以确定实际上的错误在英特尔Haswell上是否足够小，但这将是一个脆弱的优化，可能使您的代码在AMD CPU上，甚至在未来的Intel CPU上中断。 不幸的是，只需编写ISA规范就可以保证相对误差高达1.5 * 2 -12需要更多指令。 我没有看到NR迭代的任何技巧。
如果你的上限较小（如20），你可以做isqrt = static_cast 。 20*20你会得20，但20*20-1总是19。
; define test_approx_sqrt(x, off) { local s=x*x+off; local sq=s/sqrt(s); local sq_1=(s-1)/sqrt(s-1); local e=1.5*2^(-12); print sq, sq_1; print sq*e, sq_1*e; } ; test_approx_sqrt(20, 0.5) ~20.01249609618950056874 ~19.98749609130668473087 # (x+0.5)/sqrt(x+0.5) ~0.00732879495710064718 ~0.00731963968187500662 # relative error
注意val * (x +/- err) = val*x +/- val*err 。 IEEE FP mul生成的结果正确舍入到0.5ulp，因此这应该适用于FP相对错误。
无论如何，我认为你需要一个Newton-Raphson迭代。
最好的办法是在使用rsqrt执行approx_sqrt之前将0.5添加到输入中。 避开0/0 = NaN问题，并将+/-误差范围全部推到整个数字切割点的一侧（对于我们关心的范围内的数字）。
FP min / max指令具有与FP add相同的性能，并且将在关键路径上进行。 使用add而不是max也可以解决完美正方形的结果问题，可能比正确结果低几个ulp 。
编译器输出：一个不错的起点
我用clang 3.7.1和sum_int得到了相当不错的自动向sum_int结果，带有-fno-math-errno -funsafe-math-optimizations 。 -ffinite-math-only 不是必需的（但即使使用完整的-ffast-math ，clang rsqrtps在使用rsqrtps时避免使用sqrt(0) = NaN ）。
sum_fp不会自动矢量化，即使使用完整的-ffast-math 。
然而， clang的版本遇到了与你的想法相同的问题：截断rsqrt + NR的不精确结果，可能给出错误的整数。 IDK，如果这就是为什么gcc不能自动矢量化，因为它可以使用sqrtps进行大加速而不改变结果。 （至少，只要所有浮点数介于0和INT_MAX 2之间，否则转换回整数将给出INT_MIN的“无限期”结果。（符号位置位，所有其他位清零）。这是一个-ffast-math的情况-ffast-math会破坏你的程序，除非你使用-mrecip=none或其他东西。
查看godbolt上的asm输出 ：
// autovectorizes with clang, but has rounding problems. // Note the use of sqrtf, and that floorf before truncating to int is redundant. (removed because clang doesn't optimize away the roundps) int sum_int(float vals[]){ int sum = 0; for (int i = 0; i < MAX; i++) { int thisSqrt = (int) sqrtf(vals[i]); sum += std::min(thisSqrt, 0x3F); } return sum; }
要使用内部函数手动进行矢量化，我们可以查看-fno-unroll-loops的asm输出（以保持简单）。 我打算在答案中加入这个，但后来意识到它有问题。
把它放在一起：
我认为在循环中转换为int比使用addps然后addps 。 roundps是roundps上的2-uop指令（6c延迟）（SnB / IvB中为1uop）。 更糟糕的是，两个uops都需要port1，因此它们与FP add / mul竞争。 cvttps2dq是port1的1-uop指令，具有3c延迟，然后我们可以使用整数min并添加钳位和累加，因此port5可以做些事情。 使用整数向量累加器也意味着循环携带的依赖链是1个循环，因此我们不需要展开或使用多个累加器来保持多次迭代。 对于大图片来说，较小的代码总是更好（uop缓存，L1 I缓存，分支预测器）。
只要我们没有溢出32位累加器的危险，这似乎是最好的选择。 （没有对任何事情进行基准测试或甚至测试过它）。
我没有使用sqrt(x) ~= approx_recip(approx_sqrt(x))方法，因为我不知道如何进行牛顿迭代来改进它（可能它会涉及一个除法）。 而且因为复合误差较大。
这个答案的横向总和。
完整但未经测试的版本：
#include #define MAX 4096 // 2*sqrt(x) ~= 2*x*approx_rsqrt(x), with a Newton-Raphson iteration // dividing by 2 is faster in the integer domain, so we don't do it __m256 approx_2sqrt_ps256(__m256 x) { // clang / gcc usually use -3.0 and -0.5. We could do the same by using fnmsub_ps (add 3 = subtract -3), so we can share constants __m256 three = _mm256_set1_ps(3.0f); //__m256 half = _mm256_set1_ps(0.5f); // we omit the *0.5 step __m256 nr = _mm256_rsqrt_ps( x ); // initial approximation for Newton-Raphson // 1/sqrt(x) ~= nr * (3 - x*nr * nr) * 0.5 = nr*(1.5 - x*0.5*nr*nr) // sqrt(x) = x/sqrt(x) ~= (x*nr) * (3 - x*nr * nr) * 0.5 // 2*sqrt(x) ~= (x*nr) * (3 - x*nr * nr) __m256 xnr = _mm256_mul_ps( x, nr ); __m256 three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three ); // -(xnr*nr) + 3 return _mm256_mul_ps( xnr, three_minus_muls ); } // packed int32_t: correct results for inputs from 0 to well above 63*63 __m256i isqrt256_ps(__m256 x) { __m256 offset = _mm256_set1_ps(0.5f); // or subtract -0.5, to maybe share constants with compiler-generated Newton iterations. __m256 xoff = _mm256_add_ps(x, offset); // avoids 0*Inf = NaN, and rounding error before truncation __m256 approx_2sqrt_xoff = approx_2sqrt_ps256(xoff); __m256i i2sqrtx = _mm256_cvttps_epi32(approx_2sqrt_xoff); return _mm256_srli_epi32(i2sqrtx, 1); // divide by 2 with truncation // alternatively, we could mask the low bit to zero and divide by two outside the loop, but that has no advantage unless port0 turns out to be the bottleneck } __m256i isqrt256_ps_simple_exact(__m256 x) { __m256 sqrt_x = _mm256_sqrt_ps(x); __m256i isqrtx = _mm256_cvttps_epi32(sqrt_x); return isqrtx; } int hsum_epi32_avx(__m256i x256){ __m128i xhi = _mm256_extracti128_si256(x256, 1); __m128i xlo = _mm256_castsi256_si128(x256); __m128i x = _mm_add_epi32(xlo, xhi); __m128i hl = _mm_shuffle_epi32(x, _MM_SHUFFLE(1, 0, 3, 2)); hl = _mm_add_epi32(hl, x); x = _mm_shuffle_epi32(hl, _MM_SHUFFLE(2, 3, 0, 1)); hl = _mm_add_epi32(hl, x); return _mm_cvtsi128_si32(hl); } int sum_int_avx(float vals[]){ __m256i sum = _mm256_setzero_si256(); __m256i upperlimit = _mm256_set1_epi32(0x3F); for (int i = 0; i < MAX; i+=8) { __m256 v = _mm256_loadu_ps(vals+i); __m256i visqrt = isqrt256_ps(v); // assert visqrt == isqrt256_ps_simple_exact(v) or something visqrt = _mm256_min_epi32(visqrt, upperlimit); sum = _mm256_add_epi32(sum, visqrt); } return hsum_epi32_avx(sum); }
将godbolt编译成漂亮的代码 ，但我还没有测试过它。 clang使用更好的代码，gcc：clang使用4B位置的广播加载作为set1常量，而不是在编译时将它们重复为32B常量。 gcc还有一个奇怪的movdqa来复制寄存器。
无论如何，整个循环最多只有9个向量指令，而编译器生成的sum_int版本则为12。 它可能没有注意到当你将结果乘以x或类似的东西时，在Newton-Raphson迭代公式中出现的x*initial_guess(x)公共子表达式。 它还执行额外的mulps而不是psrld，因为它在转换为int之前执行* 0.5。 这就是额外的两个mulps指令的来源，并且有cmpps / blendvps。
sum_int_avx(float*): vpxor ymm3, ymm3, ymm3 xor eax, eax vbroadcastss ymm0, dword ptr [rip + .LCPI4_0] ; set1(0.5) vbroadcastss ymm1, dword ptr [rip + .LCPI4_1] ; set1(3.0) vpbroadcastd ymm2, dword ptr [rip + .LCPI4_2] ; set1(63) LBB4_1: ; latencies vaddps ymm4, ymm0, ymmword ptr [rdi + 4*rax] ; 3c vrsqrtps ymm5, ymm4 ; 7c vmulps ymm4, ymm4, ymm5 ; x*nr ; 5c vfnmadd213ps ymm5, ymm4, ymm1 ; 5c vmulps ymm4, ymm4, ymm5 ; 5c vcvttps2dq ymm4, ymm4 ; 3c vpsrld ymm4, ymm4, 1 ; 1c this would be a mulps (but not on the critical path) if we did this in the FP domain vpminsd ymm4, ymm4, ymm2 ; 1c vpaddd ymm3, ymm4, ymm3 ; 1c ; ... (those 9 insns repeated: loop unrolling) add rax, 16 cmp rax, 4096 jl .LBB4_1 ;... horizontal sum
IACA认为，没有展开，Haswell可以维持每4.15个周期一次迭代的吞吐量，在端口0和1上出现瓶颈。因此，您可以通过累积sqrt（x）* 2（使用_mm256_and_si256截断为偶数）来_mm256_and_si256 ，只在循环外除以2。
同样根据IACA，单次迭代的延迟是Haswell的38个周期。 我只得到31c，所以可能包括L1负载使用延迟等等。 无论如何，这意味着要使执行单元饱和，必须立即进行8次迭代的操作。 这是8 * ~14个未融合的域uops = 112个unfused-uops（或者更少有clang的展开）必须立即在飞行中。 Haswell的调度程序实际上只有60个条目，但ROB是192个条目 。 早期迭代的早期uops已经执行，因此它们只需要在ROB中跟踪，而不是在调度程序中跟踪。 然而，许多慢速微操都在每次迭代的开始。 尽管如此，仍然有理由希望这会使端口0和1饱和。除非L1缓存中的数据很热，否则缓存/内存带宽可能会成为瓶颈。
来自多个dep链的交叉操作也会更好。 当clang展开时，它将所有9条指令放在所有9条指令之前进行另一次迭代。 它使用了数量惊人的寄存器，因此可以将2或4次迭代的指令混合在一起。 这是编译器应该擅长的东西，但这对人类来说很麻烦。 ：/
如果编译器选择单寄存器寻址模式，那么它的效率也会稍高，因此负载可以与vaddps微融合。 gcc这样做。