assembly 如何使用AVX-512实现向量化“exp”和“log”base-2函数

p3rjfoxz  于 2023-05-23  发布在  其他
关注(0)|答案(1)|浏览(149)

对于一个游戏,我的工作-需要的能力,以运行百万的“exp”调用向量。基本上

void vector_exp(const double *x, const double *result, int n)
{
    for (int i=0 ; i<n ; i++) result[i] = exp(x[i]) ;
}

对于我的特定情况,输入都是-50..+50。需要双精度,8位小数匹配当前“exp”-以通过测试用例。
我对“日志”也有同样的挑战。输入范围为1 e-7至1 e7。
希望利用AVX 512指令-这应该能够做(理论上)8倍精度的时间。我已经检索到了glibc代码(包括为AVX构建的“C”版本和“.S”版本),但我不确定如何继续前进。
https://github.com/bminor/glibc/tree/master/sysdeps/x86_64/fpu

qyyhg6bp

qyyhg6bp1#

我相信其他的答案比我的好-运行一个非常快速和肮脏的多项式近似,我最终得到了这些。

inline __m512d exp2(const __m512d x) {
    const __m512d a = _mm512_set1_pd(0.000217549227054);
    const __m512d b = _mm512_set1_pd(0.00124218531444);
    const __m512d c = _mm512_set1_pd(0.00968102455999);
    const __m512d d = _mm512_set1_pd(0.0554821818101);
    const __m512d e = _mm512_set1_pd(0.240230073528);
    const __m512d f = _mm512_set1_pd(0.693146979806);
    const __m512d g = _mm512_set1_pd(1.0);
    const __m512d fx = _mm512_floor_pd(x);  // integer part
    const __m512d X = _mm512_sub_pd(x, fx); // fractional part
    __m512d y = _mm512_fmadd_pd(a, X, b);
    y = _mm512_fmadd_pd(y, X, c);
    y = _mm512_fmadd_pd(y, X, d);
    y = _mm512_fmadd_pd(y, X, e);
    y = _mm512_fmadd_pd(y, X, f);
    y = _mm512_fmadd_pd(y, X, g);      // polynomial approximation over [0,1)
    return _mm512_scalef_pd(y, fx);    // scale by 2^integer
}

inline __m512d exp(const __m512d x) {
    return exp2(_mm512_mul_pd(x, _mm512_set1_pd(1.442695040888963387)));
}
inline __m512d log2(const __m512d x) {
    const __m512d m = _mm512_getmant_pd(x, _MM_MANT_NORM_1_2, _MM_MANT_SIGN_zero);
    const __m512d a = _mm512_set1_pd(0.0146498917256);
    const __m512d b = _mm512_set1_pd(-0.178725976271);
    const __m512d c = _mm512_set1_pd(0.953841083567);
    const __m512d d = _mm512_set1_pd(-2.92298892586);
    const __m512d e = _mm512_set1_pd(5.68725545823);
    const __m512d f = _mm512_set1_pd(-7.4092580291);
    const __m512d g = _mm512_set1_pd(7.09194627711);
    const __m512d h = _mm512_set1_pd(-3.23671917705);
    __m512d y = _mm512_fmadd_pd(a, m, b);
    y = _mm512_fmadd_pd(y, m, c);
    y = _mm512_fmadd_pd(y, m, d);
    y = _mm512_fmadd_pd(y, m, e);
    y = _mm512_fmadd_pd(y, m, f);
    y = _mm512_fmadd_pd(y, m, g);
    y = _mm512_fmadd_pd(y, m, h);  // poly approximation over [1,2) mantissa
    return _mm512_add_pd(y, _mm512_getexp_pd(x));
}

inline __m512d log(const __m512d x) {
    return _mm512_mul_pd(log2(x), _mm512_set1_pd(0.693147180559945286));
}

跨独立的exp2()log2()操作的乱序执行可以使用六阶多项式的霍纳规则来处理FMA依赖链。
另请参阅Agner Fog的VCL实现,其目标是高精度,接近double的全精度:

  • 双精度exp_d模板,支持2.0exp(x-1)exp(x)。(有关正确的模板参数,请参阅exp2调用程序)。

使用13阶泰勒级数,代码中的注解表明它比使用Pade展开的替代版本更快:两个多项式之比。对于吞吐量来说,每多个FMA进行一次除法并不是一个灾难,特别是如果你有很多周围的代码也对每个结果进行一些非除法工作,但是这样做可能会使每个FMA进行太多的除法。

  • 双精度log_d模板。这确实使用了尾数的5阶多项式的比率。模板参数支持log(x)log(x+1)以避免丢失精度。它只做自然对数(基数为e),因此您需要将结果缩放1/ln(2)

相关问题