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

Question

下面是我想要转换的代码:该double版本的VDT的帕德精通fast_ex()约(这里的老回购资源):

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？

当然,我需要检查整体,因为我不太确定我是否已正确转换它.

编辑:我发现浮动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。 到64位整数之间的压缩转换内在函数需要AVX512DQ。

要将舍入舍入到最接近的整数，理想情况下，您应使用SSE4.1 _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC)（或将截断取零，或将floor或ceil取负-+ Inf）。

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

标量代码以结束，int n并且double px都表示相同的数值。它使用坏/笨拙的floor(val+0.5)习惯用法而不是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版本。（因为您不在乎双精度数太大而无法容纳整数的情况。）

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

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

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

或添加一个非常大的数字并再次减去它：足够大以使每个数字彼此double相距1个整数，因此常规FP舍入可以满足您的要求。（这适用于不适合32位但不double大于2 ^ 52的输入。因此无论哪种方式都可以。）另请参见如何使用SSE / AVX有效执行double / int64转换？使用这个技巧。不过，我找不到关于SO的示例。

有关：

• 指数函数使用AVX的最快的实现和指数函数的执行速度最快使用SSE与其他速度/精度的折衷，对于版本_ps（压缩单精度float）。

• 使用双精度运算的快速SSE低精度指数位于频谱的另一端，但仍适用于double。

• 在现代x86_64 CPU上，AVX / SSE运算需要花费多少个时钟周期？讨论了一些现有的库，例如SVML和Agner Fog的VCL（GPL许可）。还有glibc的libmvec。

然后当然，我需要检查整个情况，因为我不确定我是否正确转换了它。

对所有2 ^ 64 double位模式进行迭代是不切实际的，这与float只有40亿个模式不同，但也许double对尾数的低32位全为零的所有s进行迭代将是一个好的开始。即与检查循环

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对向量的简单条件设置进行向量化。在每个元素的基础上if(x) y=0;，您拥有的SIMD等效于y = (condition) ? 0 : y;，而不是。 SIMD比较产生全零/全一元素的掩码，因此您可以使用它进行混合。

例如，在这种情况下，如果您具有SSE4.1，则输入cmppd并输入blendvpd输出。或仅使用SSE2，和/或不/或混合。请参见SSE内在函数进行比较（_mm_cmpeq_ps），并且_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用全零位模式表示，因此进行比较，如果在范围内，则将生成true，然后将其与输出进行比较。（将其保留不变或将其强制为+0.0）。这要好于_mm_blendv_pd，后者在大多数Intel CPU上的成本为2微克（而AVX 128位版本在Intel上始终为2微克）。而且在AMD或Skylake上也并不差。

+-Inf用有效位= 0，指数=全1的位模式表示。（有效位数中的任何其他值都表示+ -NaN。）由于输入过大，可能仍会留下非零有效位数，因此我们不能只是将比较结果与或与最终结果进行或。我认为我们需要进行常规混合，或进行一些昂贵的混合（3 uops和向量常数）。

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

希望您的编译器将对CMP以外的所有对象实际使用ANDPS等，而不是PD，因为它短了1个字节，但相同，因为它们都是按位运算。我只是写ANDPD，所以我不必在评论中对此进行解释。

您可以通过在应用结果之前合并两个修正来缩短关键路径延迟，因此只有一个混合。但是，我认为您还需要合并比较结果。

或者由于您的上限和下限大小相同，也许您可​​以比较绝对值？（屏蔽掉initial_xdo 的符号位并执行_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);

如果您关心NaN 会发生什么情况，则可以使用cmplt代替cmpgt并重新排列initial_x。选择比较为false时，将应用修正而不是true进行修正，这意味着对于-NaN或+ NaN的输入，无序比较将导致0或+ Inf。这仍然不进行NaN传播。如果您想实现此目标，则可以将其_mm_cmpunord_pd(initial_x, initial_x)与OR或fixup。

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

您可能仍然能够使用一些这种想法的只有SSE2，也许创造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报告说vfixupimmpd，Skylake-X上只有1 uop 。