代码之家 › 专栏 › 技术社区 › markzzz

如何将vdt的pade exp fast_ex()about双版本的标量代码转换为SSE2?

exp sse2 intrinsics c++

5

markzzz · 技术社区 · 7 年前

下面是我要转换的代码: double 版本 VDT's Pade Exp fast_ex() approx (这里是 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));

1 回复 | 直到 7 年前

1

4

Peter Cordes 7 年前

是的,将两个多项式分开通常比一个大的多项式能在速度和精度之间做出更好的权衡。只要有足够的工作来隐藏 divpd 吞吐量。(最新的x86 CPU具有相当不错的fp divide吞吐量。与乘法相比仍然很糟糕,但它只有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) ,(或朝零截断,或朝-+inf地板或天花板)。

但实际上我们不需要。

标量代码以 int n 和 double px 两者都表示相同的数值。它使用 bad/buggy floor(val+0.5) idiom 而不是 rint(val) 或 nearbyint(val) 四舍五入为最接近的整数,然后将该整数转换为 double 对 int (用C++的截断语义,但这并不重要,因为 双重的 值已经是一个精确的整数。)

使用SIMD内部函数,只需将其转换为32位整数并返回即可。

__m128i n  = _mm_cvtpd_epi32( _mm_mul_pd(log2e, x) );   // round to nearest
__m128d px = _mm_cvtepi32_pd( n );

使用所需模式舍入到int,然后转换回double,相当于double->双舍入,然后像标量版本那样获取int版本。(因为你不在乎双打太大而不适合内线的情况。)

cvtsd2si和si2sd指令各自是2个uops,并将32位整数随机移动到矢量的低64位中。所以要设置64位整数,需要将这些位转换为 双重的 同样,你需要洗牌。最上面的64位 n 将是零,因此我们可以使用它创建64位整数 n 与双打并列:

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

但只要有SSE2,就有解决办法。转换为32位整数并返回是一个选项:您不关心输入太小或太大。但压缩转换 双重的 和 int Intel CPU每种方式至少需要2个Uops,因此总共需要4个Uops。但是这些UOP中只有2个需要fma单元,并且您的代码在端口5上可能不会因为所有这些乘法和加法而成为瓶颈。

或者加一个非常大的数字再减去它:足够大 双重的 间距为1个整数,所以普通的fp舍入可以满足您的需要。(这适用于不适合32位的输入,但不适用于 双重的 >2 ^ 52。所以无论哪种方法都有效。)另请看 How to efficiently perform double/int64 conversions with SSE/AVX? 使用了这个技巧。不过,我找不到这样的例子。

对于最后几行,纠正超出范围输入的输入:

这个 浮动 你引用的版本完全不在范围内。这显然是最快的方法 ,如果您的输入总是在范围内,或者您不关心超出范围的输入会发生什么。

另一种更便宜的范围检查(可能仅用于调试)是通过将打包的比较结果转换为结果,将超出范围的值转换为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);

通常,可以使用无分支比较+混合对值的简单条件设置进行矢量化。 . 而不是 if(x) y=0; ,您的simd相当于 y = (condition) ? 0 : y; ,以每个元素为基础。 SIMD比较生成一个包含所有零/所有一个元素的蒙版,这样您就可以使用它来混合。

例如,在这种情况下,如果您有SSE4.1,CMPPD输入和BLENDVPD输出。或者只用SSE2,和/不和/或混合。见 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 用全零位模式表示 ,所以做一个比较,如果在范围内,将产生真的,并且输出与之相比较。(保持不变或强制为+0.0)。这比 _mm_blendv_pd 在大多数Intel CPU上,它需要2个Uops(而在Intel上,AVX 128位版本总是需要2个Uops)。AMD和Skylake的情况也不差。

+-Inf 用有效位=0的位模式表示,指数=所有的。(有效位中的任何其他值表示+-NaN。)由于过大的输入可能仍然会留下非零的有效位,因此我们不能将比较结果和/或该值作为最终结果。我认为我们需要做一个常规的混合,或者一些昂贵的东西(3个Uops和一个向量常量)。

它为最终结果增加了2个延迟周期;andnpd和orpd都在关键路径上。cmppd和andpd不是这样的;它们可以与计算结果所做的任何操作并行运行。

希望您的编译器实际上会使用andps等,而不是pd来处理除cmp之外的所有内容,因为它比cmp短1个字节,但完全相同,因为它们都是按位操作。我写andpd只是为了不必在评论中解释。

在应用到结果之前,您可以通过组合两个修正值来缩短关键路径延迟。 所以你只有一种混合物。但是我认为你还需要把比较结果结合起来。

或者因为你的上下限是相同的,也许你可以比较绝对值?(遮住标志位 起始语言 并且做 _mm_cmplt_pd(abs_initial_x, _mm_set1_pd(details::EXP_LIMIT)) )但是你必须先把它归零或者设为+inf。

如果你有SSE4.1 第二阶段 你可以使用 起始语言 作为可能需要应用的链接地址信息的混合控件,因为 blendv 只关心混合控件的符号位(与所有位都需要匹配的and/andn/或版本不同)。

__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 如果你在乎会发生什么,重新安排 起始语言 作为一个南 .选择比较结果为“假”,则应用修正值而不是“真”,这意味着对于-nan或+nan的输入,无序的比较结果为0或+inf。这仍然不能进行NaN传播。你可以 _mm_cmpunord_pd(initial_x, initial_x) 或者说 fixup 如果你想做到这一点。

尤其是在天湖和AMD推土机/Ryzen上,SSE2 blendvpd 只有1个UOP,这应该很好。(VEX编码, vblendvpd 是2个Uops,有3个输入和一个单独的输出。)

您可能仍然能够使用一些只有SSE2的想法,或者创建 固定装置 对零进行比较,然后 _mm_and_pd 或 _mm_andnot_pd 与比较结果和+无穷大。

使用整数算术移位将符号位广播到 双重的 效率不高: psraq 不存在,只有 psraw/d . 只有64位元素大小的逻辑移位。

但你可以创造 固定装置 只有一个整数移位和掩码,以及一个位反转

__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

然后混合 固定装置 进入之内 结果 对于超出正常范围的输入。

便宜的支票 abs(initial_x) > details::EXP_LIMIT

如果exp算法已经平方化 起始语言 ,你可以和 EXP_LIMIT 平方。但事实并非如此, xx = x*x 只有在进行计算后才能创建 x .

如果您有AVX512F/VL, VFIXUPIMMPD 可能在这里很方便。它是为特殊情况输出来自“特殊”输入(如NaN和+-Inf、负、正或零)的函数而设计的,这样可以为这些情况节省比较。(例如,对于x=0的牛顿-拉斐逊倒数(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。