使用SSE最快实现自然指数函数

Question

我在寻找在SSE元素上运行的自然指数函数的近似值。即-__m128 exp（__ m128 x）。

我有一个快速但准确性似乎很低的实现：

static inline __m128 FastExpSse(__m128 x)
{
    __m128 a = _mm_set1_ps(12102203.2f); // (1 << 23) / ln(2)
    __m128i b = _mm_set1_epi32(127 * (1 << 23) - 486411);
    __m128  m87 = _mm_set1_ps(-87);
    // fast exponential function, x should be in [-87, 87]
    __m128 mask = _mm_cmpge_ps(x, m87);

    __m128i tmp = _mm_add_epi32(_mm_cvtps_epi32(_mm_mul_ps(a, x)), b);
    return _mm_and_ps(_mm_castsi128_ps(tmp), mask);
}

任何人都能有一个精度更高但速度更快（或更快）的实现吗？

如果它是用C风格写的，我会很高兴。

谢谢你。

Answer 1

回顾我那时候的笔记，我确实探索了不使用除法提高精确度的方法。我使用了相同的重新解释为浮点的技巧，但是对尾数应用了多项式校正，尾数基本上是用16位定点算法计算的(当时唯一快速的方法)。

立方分别。四分版本给你4分。精度的5位有效数字。没有必要再增加阶数，因为低精度算法的噪声会开始淹没多项式近似的误差。以下是纯C版本：

#include <stdint.h>

float fastExp3(register float x)  // cubic spline approximation
{
    union { float f; int32_t i; } reinterpreter;

    reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
    int32_t m = (reinterpreter.i >> 7) & 0xFFFF;  // copy mantissa
    // empirical values for small maximum relative error (8.34e-5):
    reinterpreter.i +=
         ((((((((1277*m) >> 14) + 14825)*m) >> 14) - 79749)*m) >> 11) - 626;
    return reinterpreter.f;
}

float fastExp4(register float x)  // quartic spline approximation
{
    union { float f; int32_t i; } reinterpreter;

    reinterpreter.i = (int32_t)(12102203.0f*x) + 127*(1 << 23);
    int32_t m = (reinterpreter.i >> 7) & 0xFFFF;  // copy mantissa
    // empirical values for small maximum relative error (1.21e-5):
    reinterpreter.i += (((((((((((3537*m) >> 16)
        + 13668)*m) >> 18) + 15817)*m) >> 14) - 80470)*m) >> 11);
    return reinterpreter.f;
}

四分位数服从（fastExp4（0f） == 1f），这对于定点迭代算法可能很重要。

这些整数乘移加序列在SSE中的效率如何？在浮点运算同样快的架构上，可以使用它来代替，减少算术噪声。这基本上会产生上面@njuffa答案的三次和四次扩展。

Answer 2

通过使用 FastExpSse（x/2）/FastExpSse（-x/2） 而不是 FastExpSse（x），可以以整数减法和浮点除法为代价，在我的算法（在上面的答案中实现 FastExpSse）中获得准确性的良好提高。这里的诀窍是将 shift 参数（298765上文）设置为零，以便分子和分母中的分段线性近似值对齐，从而为您提供实质性的误差消除。将其滚动到单个函数中：

__m128 BetterFastExpSse (__m128 x)
{
  const __m128 a = _mm_set1_ps ((1 << 22) / float(M_LN2));  // to get exp(x/2)
  const __m128i b = _mm_set1_epi32 (127 * (1 << 23));       // NB: zero shift!
  __m128i r = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
  __m128i s = _mm_add_epi32 (b, r);
  __m128i t = _mm_sub_epi32 (b, r);
  return _mm_div_ps (_mm_castsi128_ps (s), _mm_castsi128_ps (t));
}

（我不是一个硬件专家-这个部门的性能杀手有多糟糕？）

如果您需要exp（x）只是为了得到y=tanh（x）（例如，对于神经网络），请使用具有零位移位的FastExplSse，如下所示：

a = FastExpSse(x);
b = FastExpSse(-x);
y = (a - b)/(a + b);

以获得相同类型的错误消除好处。逻辑函数的工作方式类似，使用零移位的FastExPse（x/2）/（FastExPse（x/2）FastExPse（-x/2））。（这只是为了说明原理——您显然不想在这里多次评估FastExSse，而是按照上面的BetterFastExPse将其滚动到一个函数中。）

我确实从中开发了一系列高阶近似，更准确，但也更慢。未发布，但如果有人想给他们一个旋转，很乐意合作。

最后，为了好玩：在倒车档中使用FastLogSse。将其与FastOutSse链接起来，可以为您提供操作员和错误消除，并弹出一个极快的幂函数......

Answer 3

下面的C代码是我在以前回答类似问题时使用的一个算法的SSE内部函数的翻译。

基本思想是将标准指数函数的计算转换为2的幂的计算：exf（x）=ex2f（x/logf（2.0f））=ex2f（x*1.44269504）。我们将t=x*1.44269504拆分为整数i和分数f，这样t=i f和0

SSE 实现中存在的一个问题是，我们想要计算 i = floorf （t），但没有快速的方法来计算 floor（） 函数。但是，我们观察到，对于正数，floor（x） == trunc（x），对于负数，floor（x） == trunc（x） - 1，除非 x 是负整数。但是，由于核心近似可以处理 f 值 1.0f，因此对负参数使用近似值是无害的。SSE 提供了一条指令，用于将单精度浮点操作数转换为带截断的整数，因此此解决方案非常有效。

Peter Cordes指出，SSE4.1支持快速楼层功能＜code＞_mm_floor_ps（）＜/code＞，因此使用SSE4.1的变体如下所示。当启用SSE 4.1代码生成时，并非所有工具链都会自动预定义宏＜code＞__SSE4_1_。

编译器资源管理器（Godbolt）表明，gcc 7.2将下面的代码编译为16条普通SSE指令和12条SSE 4.1指令。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <emmintrin.h>
#ifdef __SSE4_1__
#include <smmintrin.h>
#endif

/* max. rel. error = 1.72863156e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 t, f, e, p, r;
    __m128i i, j;
    __m128 l2e = _mm_set1_ps (1.442695041f);  /* log2(e) */
    __m128 c0  = _mm_set1_ps (0.3371894346f);
    __m128 c1  = _mm_set1_ps (0.657636276f);
    __m128 c2  = _mm_set1_ps (1.00172476f);

    /* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */   
    t = _mm_mul_ps (x, l2e);             /* t = log2(e) * x */
#ifdef __SSE4_1__
    e = _mm_floor_ps (t);                /* floor(t) */
    i = _mm_cvtps_epi32 (e);             /* (int)floor(t) */
#else /* __SSE4_1__*/
    i = _mm_cvttps_epi32 (t);            /* i = (int)t */
    j = _mm_srli_epi32 (_mm_castps_si128 (x), 31); /* signbit(t) */
    i = _mm_sub_epi32 (i, j);            /* (int)t - signbit(t) */
    e = _mm_cvtepi32_ps (i);             /* floor(t) ~= (int)t - signbit(t) */
#endif /* __SSE4_1__*/
    f = _mm_sub_ps (t, e);               /* f = t - floor(t) */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= 2^f */
    j = _mm_slli_epi32 (i, 23);          /* i << 23 */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}

int main (void)
{
    union {
        float f[4];
        unsigned int i[4];
    } arg, res;
    double relerr, maxrelerr = 0.0;
    int i, j;
    __m128 x, y;

    float start[2] = {-0.0f, 0.0f};
    float finish[2] = {-87.33654f, 88.72283f};

    for (i = 0; i < 2; i++) {

        arg.f[0] = start[i];
        arg.i[1] = arg.i[0] + 1;
        arg.i[2] = arg.i[0] + 2;
        arg.i[3] = arg.i[0] + 3;
        do {
            memcpy (&x, &arg, sizeof(x));
            y = fast_exp_sse (x);
            memcpy (&res, &y, sizeof(y));
            for (j = 0; j < 4; j++) {
                double ref = exp ((double)arg.f[j]);
                relerr = fabs ((res.f[j] - ref) / ref);
                if (relerr > maxrelerr) {
                    printf ("arg=% 15.8e  res=%15.8e  ref=%15.8e  err=%15.8e\n", 
                            arg.f[j], res.f[j], ref, relerr);
                    maxrelerr = relerr;
                }
            }   
            arg.i[0] += 4;
            arg.i[1] += 4;
            arg.i[2] += 4;
            arg.i[3] += 4;
        } while (fabsf (arg.f[3]) < fabsf (finish[i]));
    }
    printf ("maximum relative errror = %15.8e\n", maxrelerr);
    return EXIT_SUCCESS;
}

fast_sse_exp（） 的替代设计使用众所周知的技术，在舍入到最接近的模式下提取调整后的参数 x / log（2） 的整数部分，使用众所周知的技术，将“魔术”转换常数 1.5 * 223 相加以在正确的位位置强制舍入，然后再次减去相同的数字。这要求在添加期间生效的 SSE 舍入模式是“舍入到最接近或偶数”，这是默认设置。wim在评论中指出，一些编译器在使用主动优化时，可能会将转换常数cvt的加减优化为冗余，从而干扰此代码序列的功能，因此建议检查生成的机器代码。计算 2f 的近似区间现在以零为中心，因为 -0.5

/* max. rel. error <= 1.72860465e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 t, f, p, r;
    __m128i i, j;

    const __m128 l2e = _mm_set1_ps (1.442695041f); /* log2(e) */
    const __m128 cvt = _mm_set1_ps (12582912.0f);  /* 1.5 * (1 << 23) */
    const __m128 c0 =  _mm_set1_ps (0.238428936f);
    const __m128 c1 =  _mm_set1_ps (0.703448006f);
    const __m128 c2 =  _mm_set1_ps (1.000443142f);

    /* exp(x) = 2^i * 2^f; i = rint (log2(e) * x), -0.5 <= f <= 0.5 */
    t = _mm_mul_ps (x, l2e);             /* t = log2(e) * x */
    r = _mm_sub_ps (_mm_add_ps (t, cvt), cvt); /* r = rint (t) */
    f = _mm_sub_ps (t, r);               /* f = t - rint (t) */
    i = _mm_cvtps_epi32 (t);             /* i = (int)t */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= exp2(f) */
    j = _mm_slli_epi32 (i, 23);          /* i << 23 */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}

问题中代码的算法似乎取自Nicol N.Schraudolph的工作，他巧妙地利用了IEEE-754二进制浮点格式的半对数特性：

N、 施劳多夫。“指数函数的快速紧凑逼近”，《神经计算》，第11（4）期，1999年5月，第853-862页。

删除参数钳制代码后，它减少到只有三个 SSE 指令。486411“神奇”校正常数并不是在整个输入域内最小化最大相对误差的最佳选择。基于简单的二进制搜索，值298765似乎更优越，将FastExpSse（）的最大相对误差降低到3.56e-2，而fast_exp_sse（）的最大相对误差为1.73e-3。

/* max. rel. error = 3.55959567e-2 on [-87.33654, 88.72283] */
__m128 FastExpSse (__m128 x)
{
    __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
    __m128i b = _mm_set1_epi32 (127 * (1 << 23) - 298765);
    __m128i t = _mm_add_epi32 (_mm_cvtps_epi32 (_mm_mul_ps (a, x)), b);
    return _mm_castsi128_ps (t);
}

Schraudolph的算法基本上使用线性近似2f ~= 1.0 f表示[0，1]中的f，其精度可以通过添加二次项来提高。Schraudolph方法的聪明部分是计算2i * 2f，而不明确地将整数部分i = floor（x * 1.44269504）与分数分开。我认为没有办法将这个技巧扩展到二次近似，但人们当然可以将Schraudolph的floful（）计算与上面使用的二次近似结合起来：

/* max. rel. error <= 1.72886892e-3 on [-87.33654, 88.72283] */
__m128 fast_exp_sse (__m128 x)
{
    __m128 f, p, r;
    __m128i t, j;
    const __m128 a = _mm_set1_ps (12102203.0f); /* (1 << 23) / log(2) */
    const __m128i m = _mm_set1_epi32 (0xff800000); /* mask for integer bits */
    const __m128 ttm23 = _mm_set1_ps (1.1920929e-7f); /* exp2(-23) */
    const __m128 c0 = _mm_set1_ps (0.3371894346f);
    const __m128 c1 = _mm_set1_ps (0.657636276f);
    const __m128 c2 = _mm_set1_ps (1.00172476f);

    t = _mm_cvtps_epi32 (_mm_mul_ps (a, x));
    j = _mm_and_si128 (t, m);            /* j = (int)(floor (x/log(2))) << 23 */
    t = _mm_sub_epi32 (t, j);
    f = _mm_mul_ps (ttm23, _mm_cvtepi32_ps (t)); /* f = (x/log(2)) - floor (x/log(2)) */
    p = c0;                              /* c0 */
    p = _mm_mul_ps (p, f);               /* c0 * f */
    p = _mm_add_ps (p, c1);              /* c0 * f + c1 */
    p = _mm_mul_ps (p, f);               /* (c0 * f + c1) * f */
    p = _mm_add_ps (p, c2);              /* p = (c0 * f + c1) * f + c2 ~= 2^f */
    r = _mm_castsi128_ps (_mm_add_epi32 (j, _mm_castps_si128 (p))); /* r = p * 2^i*/
    return r;
}