如何将 VDT 的 Pade Exp fast_ex() 的双重版本的标量代码大约转换为 SSE2？
How to convert scalar code of the double version of VDT's Pade Exp fast_ex() approx into SSE2?
这是我要转换的代码：VDT's Pade Exp fast_ex() approx (here's the old repo 资源的 double 版本）：
inline double fast_exp(double initial_x){
double x = initial_x;
double px=details::fpfloor(details::LOG2E * x +0.5);
const int32_t n = int32_t(px);
x -= px * 6.93145751953125E-1;
x -= px * 1.42860682030941723212E-6;
const double xx = x * x;
// px = x * P(x**2).
px = details::PX1exp;
px *= xx;
px += details::PX2exp;
px *= xx;
px += details::PX3exp;
px *= x;
// Evaluate Q(x**2).
double qx = details::QX1exp;
qx *= xx;
qx += details::QX2exp;
qx *= xx;
qx += details::QX3exp;
qx *= xx;
qx += details::QX4exp;
// e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
x = px / (qx - px);
x = 1.0 + 2.0 * x;
// Build 2^n in double.
x *= details::uint642dp(( ((uint64_t)n) +1023)<<52);
if (initial_x > details::EXP_LIMIT)
x = std::numeric_limits<double>::infinity();
if (initial_x < -details::EXP_LIMIT)
x = 0.;
return x;
}
我知道了：
__m128d PExpSSE_dbl(__m128d x) {
__m128d initial_x = x;
__m128d half = _mm_set1_pd(0.5);
__m128d one = _mm_set1_pd(1.0);
__m128d log2e = _mm_set1_pd(1.4426950408889634073599);
__m128d p1 = _mm_set1_pd(1.26177193074810590878E-4);
__m128d p2 = _mm_set1_pd(3.02994407707441961300E-2);
__m128d p3 = _mm_set1_pd(9.99999999999999999910E-1);
__m128d q1 = _mm_set1_pd(3.00198505138664455042E-6);
__m128d q2 = _mm_set1_pd(2.52448340349684104192E-3);
__m128d q3 = _mm_set1_pd(2.27265548208155028766E-1);
__m128d q4 = _mm_set1_pd(2.00000000000000000009E0);
__m128d px = _mm_add_pd(_mm_mul_pd(log2e, x), half);
__m128d t = _mm_cvtepi64_pd(_mm_cvttpd_epi64(px));
px = _mm_sub_pd(t, _mm_and_pd(_mm_cmplt_pd(px, t), one));
__m128i n = _mm_cvtpd_epi64(px);
x = _mm_sub_pd(x, _mm_mul_pd(px, _mm_set1_pd(6.93145751953125E-1)));
x = _mm_sub_pd(x, _mm_mul_pd(px, _mm_set1_pd(1.42860682030941723212E-6)));
__m128d xx = _mm_mul_pd(x, x);
px = _mm_mul_pd(xx, p1);
px = _mm_add_pd(px, p2);
px = _mm_mul_pd(px, xx);
px = _mm_add_pd(px, p3);
px = _mm_mul_pd(px, x);
__m128d qx = _mm_mul_pd(xx, q1);
qx = _mm_add_pd(qx, q2);
qx = _mm_mul_pd(xx, qx);
qx = _mm_add_pd(qx, q3);
qx = _mm_mul_pd(xx, qx);
qx = _mm_add_pd(qx, q4);
x = _mm_div_pd(px, _mm_sub_pd(qx, px));
x = _mm_add_pd(one, _mm_mul_pd(_mm_set1_pd(2.0), x));
n = _mm_add_epi64(n, _mm_set1_epi64x(1023));
n = _mm_slli_epi64(n, 52);
// return?
}
但我无法完成最后几行 - 即此代码：
if (initial_x > details::EXP_LIMIT)
x = std::numeric_limits<double>::infinity();
if (initial_x < -details::EXP_LIMIT)
x = 0.;
return x;
您将如何在 SSE2 中进行转换？
当然我需要检查整体，因为我不太确定我是否正确转换了它。
编辑：我找到了 float exp 的 SSE 转换 - 即来自：
/* multiply by power of 2 */
z *= details::uint322sp((n + 0x7f) << 23);
if (initial_x > details::MAXLOGF) z = std::numeric_limits<float>::infinity();
if (initial_x < details::MINLOGF) z = 0.f;
return z;
对此：
n = _mm_add_epi32(n, _mm_set1_epi32(0x7f));
n = _mm_slli_epi32(n, 23);
return _mm_mul_ps(z, _mm_castsi128_ps(n));
是的，与一个巨大的多项式相比，相除两个多项式通常可以在速度和精度之间取得更好的权衡。只要有足够的工作来隐藏 divpd 吞吐量。 （最新的 x86 CPU 具有相当不错的 FP 除法吞吐量。与乘法相比仍然很糟糕，但它只有 1 uop，所以如果你很少使用它，它不会停止管道，即混合大量乘法。包括在周围的代码中使用 exp)
但是，_mm_cvtepi64_pd(_mm_cvttpd_epi64(px)); 不适用于 SSE2。 Packed-conversion intrinsics to/from 64-bit integers requires AVX512DQ.
要压缩舍入到最接近的整数，理想情况下您会使用 SSE4.1 _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC)（或 t运行cation 趋向于零，或 floor 或 ceil 趋向于 -+Inf）。
但我们实际上并不需要它。
标量代码以 int n 和 double px 结束，两者都表示相同的数值。它使用 bad/buggy floor(val+0.5) idiom 而不是 rint(val) 或 nearbyint(val) 舍入到最接近的值，然后将已经是整数的 double 转换为 int （使用 C++ 的 t运行阳离子语义，但这并不重要，因为 double 值已经是一个精确的整数。）
使用 SIMD 内部函数，似乎最简单的方法就是转换为 32 位整数并返回。
__m128i n = _mm_cvtpd_epi32( _mm_mul_pd(log2e, x) ); // round to nearest
__m128d px = _mm_cvtepi32_pd( n );
使用所需模式舍入为 int，然后转换回 double，相当于 double->double 舍入，然后像标量版本一样获取它的 int 版本。 （因为你不关心双打太大而无法放入 int 会发生什么。）
cvtsd2si 和 si2sd 指令各为 2 微指令，并将 32 位整数打乱以打包到向量的低 64 位中。因此，要设置 64 位整数移位以再次将这些位填充到 double 中，您需要洗牌。 n 的前 64 位将为零，因此我们可以使用它来创建与双精度对齐的 64 位整数 n：
n = _mm_shuffle_epi32(n, _MM_SHUFFLE(3,1,2,0)); // 64-bit integers
但是只要使用 SSE2，就有变通办法。转换为 32 位整数并返回是一种选择：您不关心输入太小或太大。但是 double 和 int 之间的打包转换在 Intel CPU 上每路至少花费 2 微指令，所以总共 4 微指令。但是这些微指令中只有 2 微指令需要 FMA 单元，而您的代码可能不需要端口 5 上的所有这些乘法和加法都没有瓶颈。
或者加上一个非常大的数字然后再减去它：足够大以至于每个 double 相隔 1 个整数，所以正常的 FP 舍入可以满足您的要求。 （这适用于不适合 32 位的输入，但不适用于 double > 2^52。因此无论哪种方式都可行。）另请参阅使用该技巧的 。不过，我找不到关于 SO 的示例。
相关：
and 具有其他速度/精度权衡的版本，对于 _ps（压缩单精度 float）。
处于光谱的另一端，但仍然适用于 double.
讨论了一些现有的库，例如 SVML 和 Agner Fog 的 VCL（已获得 GPL 许可）。而glibc的libmvec.
Then of course I need to check the whole, since I'm not quite sure I've converted it correctly.
迭代所有 2^64 double 位模式是不切实际的，不像 float 只有 40 亿，但也许迭代所有 double 具有他们尾数的低 32 位全为零将是一个好的开始。即检查循环
bitpatterns = _mm_add_epi64(bitpatterns, _mm_set1_epi64x( 1ULL << 32 ));
doubles = _mm_castsi128_pd(bitpatterns);
https://randomascii.wordpress.com/2014/01/27/theres-only-four-billion-floatsso-test-them-all/
对于最后几行，更正超出范围的输入：
您引用的 float 版本完全没有进行范围检查。这显然是最快的方法，如果您的输入总是在范围内，或者如果您不关心超出范围的输入会发生什么。
另一种更便宜的范围检查（可能仅用于调试）是通过将打包比较结果与结果进行 OR 运算，将超出范围的值转换为 NaN。 （全 1 位模式表示 NaN。）
__m128d out_of_bounds = _mm_cmplt_pd( limit, abs(initial_x) ); // abs = mask off the sign bit
result = _mm_or_pd(result, out_of_bounds);
一般来说，您可以使用无分支比较+混合向量化一个值的简单条件设置。在每个元素的基础上，您拥有 y = (condition) ? 0 : y; 的 SIMD 等效项，而不是 if(x) y=0;。 SIMD 比较生成全零/全一元素的掩码，因此您可以使用它进行混合。
例如在这种情况下，如果您有 SSE4.1，则 cmppd 输入和 blendvpd 输出。或者仅使用 SSE2，and/andnot/or 进行混合。请参阅 SSE intrinsics for comparison (_mm_cmpeq_ps) and assignment operation 了解两者的 _ps 版本，_pd 是相同的。
在 asm 中它看起来像这样：
; result in xmm0 (in need of fixups for out of range inputs)
; initial_x in xmm2
; constants:
; xmm5 = limit
; xmm6 = +Inf
cmpltpd xmm2, xmm5 ; xmm2 = input_x < limit ? 0xffff... : 0
andpd xmm0, xmm2 ; result = result or 0
andnpd xmm2, xmm6 ; xmm2 = 0 or +Inf (In that order because we used ANDN)
orpd xmm0, xmm2 ; result |= 0 or +Inf
; xmm0 = (input < limit) ? result : +Inf
（在早期版本的答案中，我认为我可能正在保存一个 movaps 来复制一个寄存器，但这只是一个沼泽标准混合。它破坏了 initial_x，所以编译器需要不过，在计算 result 时复制该寄存器。）
针对这种特殊情况的优化
或在这种情况下，0.0由全零位模式表示，因此如果在范围内则进行比较，结果为真，并且以及输出结果。 （保持不变或强制其为 +0.0）。这比 _mm_blendv_pd 要好，后者在大多数 Intel CPU 上花费 2 微指令（而 AVX 128 位版本在英特尔上总是花费 2 微指令）。在 AMD 或 Skylake 上也不差。
+-Inf 由有效位数=0、指数=全一的位模式表示。 （尾数中的任何其他值表示 +-NaN。）由于输入过大可能仍会留下非零尾数，我们不能只对比较结果进行 AND 运算，然后将其与最终结果进行 OR 运算。我认为我们需要进行常规混合，或者做一些昂贵的事情（3 微指令和一个向量常数）。
它为最终结果增加了2个周期的延迟； ANDNPD 和 ORPD 都在关键路径上。 CMPPD 和 ANDPD 不是；他们可以 运行 与您计算结果的任何操作并行。
希望你的编译器实际上使用 ANDPS 等等，而不是 PD，除了 CMP 之外的所有东西，因为它短 1 个字节但相同，因为它们都是按位操作。我写 ANDPD 只是为了不必在评论中解释这一点。
您可以通过在应用到结果之前组合两个修正来缩短关键路径延迟，因此您只有一个混合。但是我认为你还需要结合比较结果。
或者你的上限和下限是同一个量级，也许你可以比较一下绝对值？ （屏蔽掉 initial_x 的符号位并执行 _mm_cmplt_pd(abs_initial_x, _mm_set1_pd(details::EXP_LIMIT))）。但是你必须弄清楚是归零还是设置为+Inf。
如果您有 _mm_blendv_pd 的 SSE4.1，您可以使用 initial_x 本身作为可能需要应用的修正的混合控件，因为 blendv 只关心符号混合控制的位（与 AND/ANDN/OR 版本不同，其中所有位都需要匹配。）
__m128d fixup = _mm_blendv_pd( _mm_setzero_pd(), _mm_set1_pd(INFINITY), initial_x ); // fixup = (initial_x signbit) ? 0 : +Inf
// see below for generating fixup with an SSE2 integer arithmetic-shift
const signbit_mask = _mm_castsi128_pd(_mm_set1_epi64x(0x7fffffffffffffff)); // ~ set1(-0.0)
__m128d abs_init_x = _mm_and_pd( initial_x, signbit_mask );
__m128d out_of_range = _mm_cmpgt_pd(abs_init_x, details::EXP_LIMIT);
// Conditionally apply the fixup to result
result = _mm_blendv_pd(result, fixup, out_of_range);
可能使用 cmplt 而不是 cmpgt 并重新排列，如果你关心 initial_x 是 NaN 会发生什么。选择比较以便 false 应用修正而不是 true 将意味着对于 -NaN 或 +NaN 的输入，无序比较会导致 0 或 +Inf。这仍然不进行 NaN 传播。如果你想做到这一点，你可以 _mm_cmpunord_pd(initial_x, initial_x) 并将其或运算为 fixup。
特别是在 Skylake 和 AMD Bulldozer/Ryzen 上，其中 SSE2 blendvpd 只有 1 uop，这应该很不错。 （VEX 编码，vblendvpd 是 2 微指令，有 3 个输入和一个单独的输出。）
您可能仍然可以仅在 SSE2 上使用这个想法的一部分，也许通过与零进行比较然后 _mm_and_pd 或 _mm_andnot_pd 与比较结果创建 fixup和 +Infinity。
使用整数算术移位将符号位广播到 double 中的每个位置效率不高：psraq 不存在，只有 psraw/d。 64 位元素大小只有逻辑移位。
但是您可以创建 fixup，只需要一个整数移位和掩码，以及一个按位反转
__m128i ix = _mm_castsi128_pd(initial_x);
__m128i ifixup = _mm_srai_epi32(ix, 11); // all 11 bits of exponent field = sign bit
ifixup = _mm_and_si128(ifixup, _mm_set1_epi64x(0x7FF0000000000000ULL) ); // clear other bits
// ix = the bit pattern for 0 (non-negative x) or +Inf (negative x)
__m128d fixup = _mm_xor_si128(ifixup, _mm_set1_epi32(-1)); // bitwise invert
然后将 fixup 混合到 result 中，以正常处理超出范围的输入。
廉价检查abs(initial_x) > details::EXP_LIMIT
如果 exp 算法已经对 initial_x 求平方，您可以与 EXP_LIMIT 平方进行比较。但事实并非如此，xx = x*x 仅在经过一些计算后发生 x.
如果您有 AVX512F/VL，VFIXUPIMMPD 在这里可能会派上用场。它专为特殊情况输出来自 "special" 输入（如 NaN 和 +-Inf、负数、正数或零）的函数而设计，从而为这些情况保存比较。 （例如，在 x=0 的 Newton-Raphson 倒数 (x) 之后。）
但是你的两个特殊情况都需要比较。还是他们？
如果您对输入求平方并减去，只需花费一个 FMA 即可执行 initial_x * initial_x - details::EXP_LIMIT * details::EXP_LIMIT 以创建对 abs(initial_x) < details::EXP_LIMIT 为负的结果，否则为非负。
Agner Fog 报告说 vfixupimmpd 在 Skylake-X 上只有 1 uop。
这是我要转换的代码：VDT's Pade Exp fast_ex() approx (here's the old repo 资源的 double 版本）：
inline double fast_exp(double initial_x){
double x = initial_x;
double px=details::fpfloor(details::LOG2E * x +0.5);
const int32_t n = int32_t(px);
x -= px * 6.93145751953125E-1;
x -= px * 1.42860682030941723212E-6;
const double xx = x * x;
// px = x * P(x**2).
px = details::PX1exp;
px *= xx;
px += details::PX2exp;
px *= xx;
px += details::PX3exp;
px *= x;
// Evaluate Q(x**2).
double qx = details::QX1exp;
qx *= xx;
qx += details::QX2exp;
qx *= xx;
qx += details::QX3exp;
qx *= xx;
qx += details::QX4exp;
// e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
x = px / (qx - px);
x = 1.0 + 2.0 * x;
// Build 2^n in double.
x *= details::uint642dp(( ((uint64_t)n) +1023)<<52);
if (initial_x > details::EXP_LIMIT)
x = std::numeric_limits<double>::infinity();
if (initial_x < -details::EXP_LIMIT)
x = 0.;
return x;
}
我知道了：
__m128d PExpSSE_dbl(__m128d x) {
__m128d initial_x = x;
__m128d half = _mm_set1_pd(0.5);
__m128d one = _mm_set1_pd(1.0);
__m128d log2e = _mm_set1_pd(1.4426950408889634073599);
__m128d p1 = _mm_set1_pd(1.26177193074810590878E-4);
__m128d p2 = _mm_set1_pd(3.02994407707441961300E-2);
__m128d p3 = _mm_set1_pd(9.99999999999999999910E-1);
__m128d q1 = _mm_set1_pd(3.00198505138664455042E-6);
__m128d q2 = _mm_set1_pd(2.52448340349684104192E-3);
__m128d q3 = _mm_set1_pd(2.27265548208155028766E-1);
__m128d q4 = _mm_set1_pd(2.00000000000000000009E0);
__m128d px = _mm_add_pd(_mm_mul_pd(log2e, x), half);
__m128d t = _mm_cvtepi64_pd(_mm_cvttpd_epi64(px));
px = _mm_sub_pd(t, _mm_and_pd(_mm_cmplt_pd(px, t), one));
__m128i n = _mm_cvtpd_epi64(px);
x = _mm_sub_pd(x, _mm_mul_pd(px, _mm_set1_pd(6.93145751953125E-1)));
x = _mm_sub_pd(x, _mm_mul_pd(px, _mm_set1_pd(1.42860682030941723212E-6)));
__m128d xx = _mm_mul_pd(x, x);
px = _mm_mul_pd(xx, p1);
px = _mm_add_pd(px, p2);
px = _mm_mul_pd(px, xx);
px = _mm_add_pd(px, p3);
px = _mm_mul_pd(px, x);
__m128d qx = _mm_mul_pd(xx, q1);
qx = _mm_add_pd(qx, q2);
qx = _mm_mul_pd(xx, qx);
qx = _mm_add_pd(qx, q3);
qx = _mm_mul_pd(xx, qx);
qx = _mm_add_pd(qx, q4);
x = _mm_div_pd(px, _mm_sub_pd(qx, px));
x = _mm_add_pd(one, _mm_mul_pd(_mm_set1_pd(2.0), x));
n = _mm_add_epi64(n, _mm_set1_epi64x(1023));
n = _mm_slli_epi64(n, 52);
// return?
}
但我无法完成最后几行 - 即此代码：
if (initial_x > details::EXP_LIMIT)
x = std::numeric_limits<double>::infinity();
if (initial_x < -details::EXP_LIMIT)
x = 0.;
return x;
您将如何在 SSE2 中进行转换？
当然我需要检查整体，因为我不太确定我是否正确转换了它。
编辑：我找到了 float exp 的 SSE 转换 - 即来自：
/* multiply by power of 2 */
z *= details::uint322sp((n + 0x7f) << 23);
if (initial_x > details::MAXLOGF) z = std::numeric_limits<float>::infinity();
if (initial_x < details::MINLOGF) z = 0.f;
return z;
对此：
n = _mm_add_epi32(n, _mm_set1_epi32(0x7f));
n = _mm_slli_epi32(n, 23);
return _mm_mul_ps(z, _mm_castsi128_ps(n));
是的，与一个巨大的多项式相比，相除两个多项式通常可以在速度和精度之间取得更好的权衡。只要有足够的工作来隐藏 divpd 吞吐量。 （最新的 x86 CPU 具有相当不错的 FP 除法吞吐量。与乘法相比仍然很糟糕，但它只有 1 uop，所以如果你很少使用它，它不会停止管道，即混合大量乘法。包括在周围的代码中使用 exp)
但是，_mm_cvtepi64_pd(_mm_cvttpd_epi64(px)); 不适用于 SSE2。 Packed-conversion intrinsics to/from 64-bit integers requires AVX512DQ.
要压缩舍入到最接近的整数，理想情况下您会使用 SSE4.1 _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC)（或 t运行cation 趋向于零，或 floor 或 ceil 趋向于 -+Inf）。
但我们实际上并不需要它。
标量代码以 int n 和 double px 结束，两者都表示相同的数值。它使用 bad/buggy floor(val+0.5) idiom 而不是 rint(val) 或 nearbyint(val) 舍入到最接近的值，然后将已经是整数的 double 转换为 int （使用 C++ 的 t运行阳离子语义，但这并不重要，因为 double 值已经是一个精确的整数。）
使用 SIMD 内部函数，似乎最简单的方法就是转换为 32 位整数并返回。
__m128i n = _mm_cvtpd_epi32( _mm_mul_pd(log2e, x) ); // round to nearest
__m128d px = _mm_cvtepi32_pd( n );
使用所需模式舍入为 int，然后转换回 double，相当于 double->double 舍入，然后像标量版本一样获取它的 int 版本。 （因为你不关心双打太大而无法放入 int 会发生什么。）
cvtsd2si 和 si2sd 指令各为 2 微指令，并将 32 位整数打乱以打包到向量的低 64 位中。因此，要设置 64 位整数移位以再次将这些位填充到 double 中，您需要洗牌。 n 的前 64 位将为零，因此我们可以使用它来创建与双精度对齐的 64 位整数 n：
n = _mm_shuffle_epi32(n, _MM_SHUFFLE(3,1,2,0)); // 64-bit integers
但是只要使用 SSE2，就有变通办法。转换为 32 位整数并返回是一种选择：您不关心输入太小或太大。但是 double 和 int 之间的打包转换在 Intel CPU 上每路至少花费 2 微指令，所以总共 4 微指令。但是这些微指令中只有 2 微指令需要 FMA 单元，而您的代码可能不需要端口 5 上的所有这些乘法和加法都没有瓶颈。
或者加上一个非常大的数字然后再减去它：足够大以至于每个 double 相隔 1 个整数，所以正常的 FP 舍入可以满足您的要求。 （这适用于不适合 32 位的输入，但不适用于 double > 2^52。因此无论哪种方式都可行。）另请参阅使用该技巧的
相关：
_ps（压缩单精度float）。double.libmvec.
Then of course I need to check the whole, since I'm not quite sure I've converted it correctly.
迭代所有 2^64 double 位模式是不切实际的，不像 float 只有 40 亿，但也许迭代所有 double 具有他们尾数的低 32 位全为零将是一个好的开始。即检查循环
bitpatterns = _mm_add_epi64(bitpatterns, _mm_set1_epi64x( 1ULL << 32 ));
doubles = _mm_castsi128_pd(bitpatterns);
https://randomascii.wordpress.com/2014/01/27/theres-only-four-billion-floatsso-test-them-all/
对于最后几行，更正超出范围的输入：
您引用的 float 版本完全没有进行范围检查。这显然是最快的方法，如果您的输入总是在范围内，或者如果您不关心超出范围的输入会发生什么。
另一种更便宜的范围检查（可能仅用于调试）是通过将打包比较结果与结果进行 OR 运算，将超出范围的值转换为 NaN。 （全 1 位模式表示 NaN。）
__m128d out_of_bounds = _mm_cmplt_pd( limit, abs(initial_x) ); // abs = mask off the sign bit
result = _mm_or_pd(result, out_of_bounds);
一般来说，您可以使用无分支比较+混合向量化一个值的简单条件设置。在每个元素的基础上，您拥有 y = (condition) ? 0 : y; 的 SIMD 等效项，而不是 if(x) y=0;。 SIMD 比较生成全零/全一元素的掩码，因此您可以使用它进行混合。
例如在这种情况下，如果您有 SSE4.1，则 cmppd 输入和 blendvpd 输出。或者仅使用 SSE2，and/andnot/or 进行混合。请参阅 SSE intrinsics for comparison (_mm_cmpeq_ps) and assignment operation 了解两者的 _ps 版本，_pd 是相同的。
在 asm 中它看起来像这样：
; result in xmm0 (in need of fixups for out of range inputs)
; initial_x in xmm2
; constants:
; xmm5 = limit
; xmm6 = +Inf
cmpltpd xmm2, xmm5 ; xmm2 = input_x < limit ? 0xffff... : 0
andpd xmm0, xmm2 ; result = result or 0
andnpd xmm2, xmm6 ; xmm2 = 0 or +Inf (In that order because we used ANDN)
orpd xmm0, xmm2 ; result |= 0 or +Inf
; xmm0 = (input < limit) ? result : +Inf
（在早期版本的答案中，我认为我可能正在保存一个 movaps 来复制一个寄存器，但这只是一个沼泽标准混合。它破坏了 initial_x，所以编译器需要不过，在计算 result 时复制该寄存器。）
针对这种特殊情况的优化
或在这种情况下，0.0由全零位模式表示，因此如果在范围内则进行比较，结果为真，并且以及输出结果。 （保持不变或强制其为 +0.0）。这比 _mm_blendv_pd 要好，后者在大多数 Intel CPU 上花费 2 微指令（而 AVX 128 位版本在英特尔上总是花费 2 微指令）。在 AMD 或 Skylake 上也不差。
+-Inf 由有效位数=0、指数=全一的位模式表示。 （尾数中的任何其他值表示 +-NaN。）由于输入过大可能仍会留下非零尾数，我们不能只对比较结果进行 AND 运算，然后将其与最终结果进行 OR 运算。我认为我们需要进行常规混合，或者做一些昂贵的事情（3 微指令和一个向量常数）。
它为最终结果增加了2个周期的延迟； ANDNPD 和 ORPD 都在关键路径上。 CMPPD 和 ANDPD 不是；他们可以 运行 与您计算结果的任何操作并行。
希望你的编译器实际上使用 ANDPS 等等，而不是 PD，除了 CMP 之外的所有东西，因为它短 1 个字节但相同，因为它们都是按位操作。我写 ANDPD 只是为了不必在评论中解释这一点。
您可以通过在应用到结果之前组合两个修正来缩短关键路径延迟，因此您只有一个混合。但是我认为你还需要结合比较结果。
或者你的上限和下限是同一个量级，也许你可以比较一下绝对值？ （屏蔽掉 initial_x 的符号位并执行 _mm_cmplt_pd(abs_initial_x, _mm_set1_pd(details::EXP_LIMIT))）。但是你必须弄清楚是归零还是设置为+Inf。
如果您有 _mm_blendv_pd 的 SSE4.1，您可以使用 initial_x 本身作为可能需要应用的修正的混合控件，因为 blendv 只关心符号混合控制的位（与 AND/ANDN/OR 版本不同，其中所有位都需要匹配。）
__m128d fixup = _mm_blendv_pd( _mm_setzero_pd(), _mm_set1_pd(INFINITY), initial_x ); // fixup = (initial_x signbit) ? 0 : +Inf
// see below for generating fixup with an SSE2 integer arithmetic-shift
const signbit_mask = _mm_castsi128_pd(_mm_set1_epi64x(0x7fffffffffffffff)); // ~ set1(-0.0)
__m128d abs_init_x = _mm_and_pd( initial_x, signbit_mask );
__m128d out_of_range = _mm_cmpgt_pd(abs_init_x, details::EXP_LIMIT);
// Conditionally apply the fixup to result
result = _mm_blendv_pd(result, fixup, out_of_range);
可能使用 cmplt 而不是 cmpgt 并重新排列，如果你关心 initial_x 是 NaN 会发生什么。选择比较以便 false 应用修正而不是 true 将意味着对于 -NaN 或 +NaN 的输入，无序比较会导致 0 或 +Inf。这仍然不进行 NaN 传播。如果你想做到这一点，你可以 _mm_cmpunord_pd(initial_x, initial_x) 并将其或运算为 fixup。
特别是在 Skylake 和 AMD Bulldozer/Ryzen 上，其中 SSE2 blendvpd 只有 1 uop，这应该很不错。 （VEX 编码，vblendvpd 是 2 微指令，有 3 个输入和一个单独的输出。）
您可能仍然可以仅在 SSE2 上使用这个想法的一部分，也许通过与零进行比较然后 _mm_and_pd 或 _mm_andnot_pd 与比较结果创建 fixup和 +Infinity。
使用整数算术移位将符号位广播到 double 中的每个位置效率不高：psraq 不存在，只有 psraw/d。 64 位元素大小只有逻辑移位。
但是您可以创建 fixup，只需要一个整数移位和掩码，以及一个按位反转
__m128i ix = _mm_castsi128_pd(initial_x);
__m128i ifixup = _mm_srai_epi32(ix, 11); // all 11 bits of exponent field = sign bit
ifixup = _mm_and_si128(ifixup, _mm_set1_epi64x(0x7FF0000000000000ULL) ); // clear other bits
// ix = the bit pattern for 0 (non-negative x) or +Inf (negative x)
__m128d fixup = _mm_xor_si128(ifixup, _mm_set1_epi32(-1)); // bitwise invert
然后将 fixup 混合到 result 中，以正常处理超出范围的输入。
廉价检查abs(initial_x) > details::EXP_LIMIT
如果 exp 算法已经对 initial_x 求平方，您可以与 EXP_LIMIT 平方进行比较。但事实并非如此，xx = x*x 仅在经过一些计算后发生 x.
如果您有 AVX512F/VL，VFIXUPIMMPD 在这里可能会派上用场。它专为特殊情况输出来自 "special" 输入（如 NaN 和 +-Inf、负数、正数或零）的函数而设计，从而为这些情况保存比较。 （例如，在 x=0 的 Newton-Raphson 倒数 (x) 之后。）
但是你的两个特殊情况都需要比较。还是他们？
如果您对输入求平方并减去，只需花费一个 FMA 即可执行 initial_x * initial_x - details::EXP_LIMIT * details::EXP_LIMIT 以创建对 abs(initial_x) < details::EXP_LIMIT 为负的结果，否则为非负。
Agner Fog 报告说 vfixupimmpd 在 Skylake-X 上只有 1 uop。