时间戳

首页 > 解决方案 > 如何将 VDT 的 Pade Exp fast_ex() 大约双版本的标量代码转换为 SSE2？

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

问题描述

这是我要转换的代码：VDT 的 Pade Exp fast_ex() 近似double版本（这是旧的 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));

标签： c++sseintrinsicssse2exp

解决方案

是的，除以两个多项式通常可以比一个巨大的多项式更好地权衡速度和精度。只要有足够的工作来隐藏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都代表相同的数值结束。它使用bad/buggyfloor(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 版本。（因为您不在乎双打太大而无法放入 int 会发生什么。）

cvtsd2si 和 si2sd 指令各为 2 个微指令，并将 32 位整数打乱以打包在向量的低 64 位中。因此，要设置 64 位整数移位以double再次将这些位填充到 a 中，您需要进行洗牌。的前 64 位n将为零，因此我们可以使用它来创建n与双精度对齐的 64 位整数：

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

但是只有 SSE2，有一些解决方法。转换为 32 位整数并返回是一种选择：您不关心输入太小或太大。但是在 Intel CPU 上进行打包转换至少需要 2 微指令，因此总共需要 4。但是这些微指令中只有 2 个需要 FMA 单元，并且您的代码在端口 5 上可能不会因为所有这些乘法而成为瓶颈，double并且int补充道。

或者添加一个非常大的数字并再次减去它：足够大，每个double相隔 1 个整数，所以正常的 FP 舍入可以满足您的需求。（这适用于不适合 32 位但不适合double> 2^52 的输入。所以无论哪种方式都可以。）另请参阅如何使用 SSE/AVX 有效地执行双精度/整数 64 转换？它使用了那个技巧。不过，我找不到关于 SO 的示例。

有关的：

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

迭代所有 2^64double位模式是不切实际的，这与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您引用的版本完全省略了范围检查。如果您的输入始终在范围内，或者您不关心超出范围的输入会发生什么，这显然是最快的方法。

另一种更便宜的范围检查（可能仅用于调试）是通过将打包比较结果与结果进行 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);

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

例如，在这种情况下，如果您有 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由全零位模式表示，因此如果在范围内，则进行比较将产生真值，并将输出与该值相加。（保持不变或强制为+0.0）。这比_mm_blendv_pd在大多数 Intel CPU 上花费 2 uops 更好（而 AVX 128 位版本在 Intel 上总是花费 2 uops）。AMD 或 Skylake 的情况也不差。

+-Inf由significand=0，exponent=all-ones的位模式表示。（有效数字中的任何其他值都表示 +-NaN。）由于太大的输入可能仍会留下非零有效数字，因此我们不能只将比较结果和 OR 合并到最终结果中。我认为我们需要做一个常规的混合，或者一些昂贵的东西（3 uops 和一个矢量常数）。

它为最终结果增加了 2 个延迟周期；ANDNPD 和 ORPD 都在关键路径上。CMPPD 和 ANDPD 不是；它们可以与您为计算结果所做的任何事情并行运行。

希望您的编译器实际上对除 CMP 之外的所有内容都使用 ANDPS 等，而不是 PD，因为它短了 1 个字节但相同，因为它们都是按位操作。我写 ANDPD 只是为了不必在评论中解释这一点。

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

或者由于您的上下限是相同的数量级，也许您可​​以比较绝对值？（屏蔽掉initial_xand do的符号位_mm_cmplt_pd(abs_initial_x, _mm_set1_pd(details::EXP_LIMIT))）。但是你必须弄清楚是设置为零还是设置为+Inf。

如果您有 SSE4.1 _mm_blendv_pd，则可以将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代替cmpgtinitial_x并重新排列。选择 compare so false 将应用修正而不是 true 将意味着对于 -NaN 或 +NaN 的输入，无序比较结果为 0 或 +Inf。这仍然不会进行 NaN 传播。如果你想让它发生，你可以_mm_cmpunord_pd(initial_x, initial_x)和 OR 成。fixup

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

您可能仍然可以仅对 SSE2 使用其中的一些想法，可能fixup通过与零进行比较然后_mm_and_pd或_mm_andnot_pd使用比较结果和 +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在这里可能会很方便。它专为特殊情况输出来自“特殊”输入（如 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 报告说vfixupimmpdSkylake-X 上只有 1 uop。

推荐阅读