下面的C代码是我在
previous answer中用于类似问题的算法的SSE内在函数的转换.
基本思想是将标准指数函数的计算转换为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
#include
#include
#include
#include
#ifdef __SSE4_1__
#include
#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