我正在寻找在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风格编写，我会很高兴。
谢谢。
下面的C代码是我在previous answer中使用的类似问题的算法的SSE内在函数的翻译。
基本思想是将标准指数函数的计算转换为2的幂的计算：expf (x) = exp2f (x / logf (2.0f)) = exp2f (x * 1.44269504)。我们将t = x * 1.44269504分为整数i和分数f，这样t = i + f和0 <= f <= 1。现在，我们可以使用多项式逼近来计算2 f，然后通过将i添加到单精度浮点结果的指数字段中，将结果缩放2 i。
SSE实现中存在的一个问题是我们要计算i = floorf (t)，但是没有快速的方法来计算floor()函数。但是，我们观察到，对于正数floor(x) == trunc(x)，对于负数floor(x) == trunc(x) - 1，除非x是负整数。但是，由于核心近似值可以处理f的1.0f值，因此对负参数使用近似值是无害的。 SSE提供了一条指令，可以将单精度浮点操作数转换为带有截断的整数，因此这种解决方案非常有效。
Peter Cordes指出SSE4.1支持快速下限功能_mm_floor_ps()，因此，下面也显示了使用SSE4.1的变体。启用SSE 4.1代码生成时，并非所有工具链都会自动预定义宏__SSE4_1__，但gcc会自动预定义宏。
[Compiler Explorer（Godbolt）显示gcc 7.2将以下代码编译为普通SSE的sixteen instructions和SSE 4.1的twelve instructions。
#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()的另一种设计，使用众所周知的添加“魔术”转换常数1.5 * 2 23]的技术，在最近舍入模式下提取调整后参数x / log(2)的整数部分。强制四舍五入到正确的位位置，然后再次减去相同的数字。这要求加法期间有效的SSE舍入模式为“舍入到最接近或什至是”，这是默认设置。 wim在注释中指出，当使用积极优化时，某些编译器可能会优化转换常数cvt的加法和减法，这会干扰此代码序列的功能，因此建议检查机器代码产生。由于-0.5 <= f <= 0.5，需要不同的核心近似，因此计算2 f的近似间隔现在以零为中心。/* 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. N. Schraudolph. "A fast, compact approximation of the exponential function." 神经计算，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算法基本上对[0,1]中的
〜=1.0 + f使用线性近似2 ff，并且可以通过添加二次项来提高其准确性。 Schraudolph方法的巧妙部分是计算2 [i] * 2 [f]，而没有明确地从分数中分离出整数部分i = floor(x * 1.44269504)。我看不到将这种技巧扩展到二次近似的方法，但是可以肯定的是可以将Schraudolph的floor()计算与上面使用的二次近似结合起来：
/* 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;
}
我的算法（在上面的答案中实现FastExpSse）的准确性有了很好的提高，而可以通过使用FastExpSse（x / 2）/ FastExpSse（-x / 2）进行整数减法和浮点除法来获得FastExpSse（x）的值。这里的技巧是将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）（例如，对于神经网络而言，请使用FastExpSse，零移位如下：
a = FastExpSse(x); b = FastExpSse(-x); y = (a - b)/(a + b);以获得相同类型的错误消除收益。逻辑函数的工作原理类似，使用具有零移位的FastExpSse（x / 2）/（FastExpSse（x / 2）+ FastExpSse（-x / 2））。 （这只是为了显示原理-您显然不想在这里多次评估FastExpSse，而是按照上面的BetterFastExpSse的方式将其滚动到单个函数中。）
我确实从中得出了一系列高阶近似值，但更精确但也更慢。未公开，但很高兴与他人合作，如果有人想旋转一下。最后，为了获得一些乐趣：使用倒档获得FastLogSse。与FastExpSse链接可以为您提供操作员和错误消除功能，并且弹出会弹出非常快速的强大功能...
#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答案的三次和四次扩展。
其tanh方程6，在第9页上与@NicSchraudolph答案非常相似
__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));
}
（（我不是硬件专家，这里的性能杀手有多糟？）
#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;
}
其tanh方程6，在第9页上与@NicSchraudolph答案非常相似