Roy*_*oyi 13 c optimization sse simd vectorization
我正在寻找在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风格写的话,我会很高兴的.
谢谢.
nju*_*ffa 13
下面的C代码是我在先前对类似问题的答案中使用的算法的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,然后通过加上单精度浮点结果的指数字段将结果缩放2 ii.
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的变体也如下所示.__SSE4_1__当启用SSE 4.1代码生成时,并非所有工具链都会自动预定义宏,但gcc会这样做.
编译器资源管理器(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*2 23迫使在正确的比特位置舍入,然后减去再次使用相同的数字.这要求在添加期间有效的SSE舍入模式是"舍入到最接近或甚至",这是默认值.wim在评论中指出,cvt当使用积极优化时,某些编译器可能会优化转换常量的加法和减法,从而干扰此代码序列的功能,因此建议检查生成的机器代码.为2的计算的近似间隔˚F现在围绕着零,因为-0.5 <= f <= 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二进制浮点格式的半对数性质:
NN Schraudolph."指数函数的快速,紧凑近似." Neural Computation,11(4),1999年5月,pp.853-862.
删除参数钳位代码后,它只减少到三个SSE指令."神奇"校正常数486411对于最小化整个输入域上的最大相对误差不是最佳的.基于简单的二进制搜索,该值298765似乎更优越,将最大相对误差降低FastExpSse()到3.56e-2,最大相对误差为1.73e-3 fast_exp_sse().
/* 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]中的线性逼近2 f ~ = 1.0 + ffor f,并且通过添加二次项可以提高其精度.Schraudolph方法的聪明部分是计算2 i*2 f而没有明确地将整数部分i = floor(x * 1.44269504)与分数分开.我认为无法将该技巧扩展到二次近似,但是当然可以将floor()Schraudolph 的计算与上面使用的二次近似相结合:
/* 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(x/2)/ FastExpSse(-x/2)代替FastExpSse,可以通过整数减法和浮点除法来获得算法精度的良好提高(在上面的答案中实现FastExpSse) (X).这里的技巧是将移位参数(上面的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进行链接可以为您提供操作员和错误消除功能,并且可以弹出超快速的电源功能......
回顾我当时的笔记，我确实探索了在不使用除法的情况下提高准确性的方法。我使用了相同的 reinterpret-as-float 技巧，但对尾数应用了多项式校正，尾数基本上是用 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 答案的三次和四次扩展。
| 归档时间： |
|
| 查看次数： |
4246 次 |
| 最近记录： |