我有一个快速但实际上准确度非常低的实现：
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风格写的话,我会很高兴的.
谢谢.
解决方法
基本思想是将标准指数函数的计算转换为2的幂的计算：expf(x)= exp2f(x / logf(2.0f))= exp2f(x * 1.44269504).我们将t = x * 1.44269504分成整数i和分数f,使得t = if和0 <= f <= 1.我们现在可以用多项式近似计算2f,然后通过添加将结果缩放2i 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会这样做.
编译器资源管理器(Godbolt)显示gcc 7.2将下面的代码编译为sixteen instructions,用于普通SSE,twelve instructions用于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()的另一种设计使用众所周知的技术添加“魔法”转换常量1.5 * 223来强制舍入,从而以舍入到最接近的模式提取调整后的参数x / log(2)的整数部分.正确的位位置,然后再次减去相同的数字.这要求在添加期间有效的SSE舍入模式是“舍入到最接近或甚至”,这是默认值. wim在评论中指出,当使用积极优化时,某些编译器可能会优化转换常量cvt的加法和减法,从而干扰此代码序列的功能,因此建议检查生成的机器代码.用于计算2f的近似间隔现在以零为中心,因为-0.5 <= f <= 0.5,需要不同的核近似.
/* max. rel. error <= 1.72860465e-3 on [-87.33654,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,c2); /* p = (c0 * f + c1) * f + c2 ~= exp2(f) */
j = _mm_slli_epi32 (i,_mm_castps_si128 (p))); /* r = p * 2^i*/
return r;
}
问题中代码的算法似乎取自Nicol N. Schraudolph的工作,它巧妙地利用了IEEE-754二进制浮点格式的半对数性质：
N. N. Schraudolph. “指数函数的快速,紧凑近似.” Neural Computation,11(4),1999年5月,pp.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,b);
return _mm_castsi128_ps (t);
}
Schraudolph算法基本上对[0,1]中的f使用线性逼近2f~ = 1.0 f,并且通过添加二次项可以提高其精度. Schraudolph方法的聪明部分是计算2i * 2f而没有明确地将整数部分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,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,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;
}
