c++ - 如何将VDT的Pade Exp fast_ex()双重版本的标量代码转换为SSE2？

Question

这是我要转换的代码:VDT's Pade Exp fast_ex() approx的double版本(这是old repo资源):

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));

Answer 1

是的，将两个多项式相除通常可以比一个庞大的多项式在速度和精度之间取得更好的折衷。只要有足够的工作来隐藏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) (或将截断符向零，或将floor或ceil移至-+ Inf)。

但是我们实际上并不需要。

标量代码以int n和double px结尾，它们都表示相同的数值。它使用bad/buggy floor(val+0.5) idiom而不是rint(val)或nearbyint(val)四舍五入到最接近的值，然后将已经整数的double转换为int(具有C++的截断语义，但这并不重要，因为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版本。 (因为您不在乎 double 数太大而无法容纳int的情况。)

cvtsd2si和si2sd指令各为2 oups，将32位整数进行混洗以打包到 vector 的低64位中。因此，要设置64位整数移位以将这些位再次填充到 double中，您需要进行洗牌。 n的高64位将为零，因此我们可以使用它来创建64位整数 n，并与 double 字对齐:

n = _mm_shuffle_epi32(n, _MM_SHUFFLE(3,1,2,0));   // 64-bit integers

但是只有SSE2，有一些解决方法。转换为32位整数然后返回是一种选择:您不必担心输入太小或太大。但是在Intel CPU上 double和 int之间的压缩转换至少要花费2 uops，因此总共要花费4 uop。但是这些uop中只有2 uops需要FMA单元，并且您的代码可能不会在端口5上成为所有这些瓶颈相乘并相加。

或添加一个非常大的数字并再次减去它:足够大以至于每个 double相隔1个整数，因此常规FP舍入可以满足您的要求。 (这适用于不适合32位的输入，但不适用于 double> 2 ^ 52。因此无论哪种方式都可以。)另请参见 How to efficiently perform double/int64 conversions with SSE/AVX?使用该技巧。不过，我找不到关于SO的示例。

有关的:

Fastest Implementation of Exponential Function Using AVX和Fastest Implementation of Exponential Function Using SSE具有其他速度/精度折衷的版本，适用于_ps(打包的单精度float)。

Fast SSE low precision exponential using double precision operations在频谱的另一端，但仍然是double。

How many clock cycles does cost AVX/SSE exponentiation on modern x86_64 CPU?讨论了一些现有的库，例如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亿个不同，但是也许迭代所有尾数的低32位都为零的 double会是一个好的开始。即与

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版本完全没有进行范围检查。如果您的输入将始终在范围内，或者您不关心超出范围的输入会发生什么，那么这显然是最快的方法。

另一种更便宜的范围检查(可能仅用于调试)是通过将打包比较结果与结果进行或运算，将超出范围的值转换为NaN。 (全位模式表示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);

通常，您可以使用无分支比较+ blend 将值的简单条件设置向量化。在每个元素的基础上，您拥有的SIMt等效于 if(x) y=0;，而不是 y = (condition) ? 0 : y;。 SIMD比较会产生全零/全一元素的掩码，因此您可以使用它进行混合。

例如在这种情况下，如果您具有SSE4.1，则输入cmppd并输入blendvpd输出。或仅使用SSE2，和/或不/或进行混合。两者的 _ps版本请参见 SSE intrinsics for comparison (_mm_cmpeq_ps) and assignment operation， _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由全零位模式表示，因此请进行比较，如果在范围内，则将生成true，并将其与输出进行比较。 (将其保留不变或将其强制为+0.0)。这比 _mm_blendv_pd更好，后者在大多数Intel CPU上的成本为2 uop(而AVX 128位版本在Intel上的价格始终为2 uop)。而且在AMD或Skylake上也并不差。
+-Inf用有效位= 0，指数=全1的位模式表示。 (有效位数中的任何其他值都表示+ -NaN。)由于输入过大可能仍会留下非零有效位数，因此我们不能只是将比较结果与或与最终结果进行“或”运算。我认为我们需要进行常规混合，或进行一些昂贵的混合(3 uops和 vector 常量)。

它为最终结果增加了2个周期的延迟； ANDNPD和ORPD都处于关键路径上。 CMPPD和ANDPD不是；它们可以与您执行的任何并行操作来计算结果。

希望您的编译器将对CMP以外的所有对象实际使用ANDPS等，而不是PD，因为它短了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是NaN 会发生什么情况，则可能使用 cmpgt而不是 initial_x并重新排列。选择比较为false时，将应用修正而不是true，这意味着对于-NaN或+ NaN的输入，无序比较将导致0或+ Inf。这仍然不进行NaN传播。如果要实现此目的，可以将 _mm_cmpunord_pd(initial_x, initial_x)和OR或为 fixup。

特别是在Skylake和AMD Bulldozer/Ryzen上，其中SSE2 blendvpd只有1 uop，这应该非常不错。 (VEX编码 vblendvpd为2 uops，具有3个输入和一个单独的输出。)

您也许仍然只能在SSE2上使用某些想法，也许通过对零进行比较，然后使用比较结果和+ Infinity进行 fixup或 _mm_and_pd来创建 _mm_andnot_pd。

使用整数算术移位将符号位广播到 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 可能会很方便。它设计用于特殊情况输出来自“特殊”输入(例如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报告说，在Skylake-X上， vfixupimmpd仅为1 uop。