你如何用 SSE2 处理 exp()？

Question

我正在编写一个代码，基本上利用 SSE2 来优化此代码：

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];

for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
    pC[sampleIndex] = exp((mMin + std::clamp(pA[sampleIndex] + pB[sampleIndex], 0.0, 1.0) * mRange) * ln2per12);
}

在此：

double *pA = a;
double *pB = b[voiceIndex];
double *pC = c[voiceIndex];

// SSE2
__m128d bound_lower = _mm_set1_pd(0.0);
__m128d bound_upper = _mm_set1_pd(1.0);
__m128d rangeLn2per12 = _mm_set1_pd(mRange * ln2per12);
__m128d minLn2per12 = _mm_set1_pd(mMin * ln2per12);

__m128d loaded_a = _mm_load_pd(pA);
__m128d loaded_b = _mm_load_pd(pB);
__m128d result = _mm_add_pd(loaded_a, loaded_b);
result = _mm_max_pd(bound_lower, result);
result = _mm_min_pd(bound_upper, result);
result = _mm_mul_pd(rangeLn2per12, result);
result = _mm_add_pd(minLn2per12, result);

double *pCEnd = pC + roundintup8(blockSize);
for (; pC < pCEnd; pA += 8, pB += 8, pC += 8) {
    _mm_store_pd(pC, result);

    loaded_a = _mm_load_pd(pA + 2);
    loaded_b = _mm_load_pd(pB + 2);
    result = _mm_add_pd(loaded_a, loaded_b);
    result = _mm_max_pd(bound_lower, result);
    result = _mm_min_pd(bound_upper, result);
    result = _mm_mul_pd(rangeLn2per12, result);
    result = _mm_add_pd(minLn2per12, result);
    _mm_store_pd(pC + 2, result);

    loaded_a = _mm_load_pd(pA + 4);
    loaded_b = _mm_load_pd(pB + 4);
    result = _mm_add_pd(loaded_a, loaded_b);
    result = _mm_max_pd(bound_lower, result);
    result = _mm_min_pd(bound_upper, result);
    result = _mm_mul_pd(rangeLn2per12, result);
    result = _mm_add_pd(minLn2per12, result);
    _mm_store_pd(pC + 4, result);

    loaded_a = _mm_load_pd(pA + 6);
    loaded_b = _mm_load_pd(pB + 6);
    result = _mm_add_pd(loaded_a, loaded_b);
    result = _mm_max_pd(bound_lower, result);
    result = _mm_min_pd(bound_upper, result);
    result = _mm_mul_pd(rangeLn2per12, result);
    result = _mm_add_pd(minLn2per12, result);
    _mm_store_pd(pC + 6, result);

    loaded_a = _mm_load_pd(pA + 8);
    loaded_b = _mm_load_pd(pB + 8);
    result = _mm_add_pd(loaded_a, loaded_b);
    result = _mm_max_pd(bound_lower, result);
    result = _mm_min_pd(bound_upper, result);
    result = _mm_mul_pd(rangeLn2per12, result);
    result = _mm_add_pd(minLn2per12, result);
}

而且我会说效果很好。但是，找不到 SSE2 的任何 exp 函数来完成操作链。

正在阅读 this，看来我需要从图书馆调用标准 exp()？

真的吗？这不是惩罚吗？还有其他办法吗？不同的内置函数？

我正在 MSVC、/arch:SSE2、/O2，生成 32 位代码。

Answer 1

最简单的方法是使用指数近似。基于此限制的一种可能情况

对于n = 256 = 2^8：

__m128d fastExp1(__m128d x)
{
   __m128d ret = _mm_mul_pd(_mm_set1_pd(1.0 / 256), x);
   ret = _mm_add_pd(_mm_set1_pd(1.0), ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   ret = _mm_mul_pd(ret, ret);
   return ret;
}

另一个想法是多项式展开。特别是泰勒级数展开：

__m128d fastExp2(__m128d x)
{
   const __m128d a0 = _mm_set1_pd(1.0);
   const __m128d a1 = _mm_set1_pd(1.0);
   const __m128d a2 = _mm_set1_pd(1.0 / 2);
   const __m128d a3 = _mm_set1_pd(1.0 / 2 / 3);
   const __m128d a4 = _mm_set1_pd(1.0 / 2 / 3 / 4);
   const __m128d a5 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5);
   const __m128d a6 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5 / 6);
   const __m128d a7 = _mm_set1_pd(1.0 / 2 / 3 / 4 / 5 / 6 / 7);

   __m128d ret = _mm_fmadd_pd(a7, x, a6);
   ret = _mm_fmadd_pd(ret, x, a5); 
   // If fma extention is not present use
   // ret = _mm_add_pd(_mm_mul_pd(ret, x), a5);
   ret = _mm_fmadd_pd(ret, x, a4);
   ret = _mm_fmadd_pd(ret, x, a3);
   ret = _mm_fmadd_pd(ret, x, a2);
   ret = _mm_fmadd_pd(ret, x, a1);
   ret = _mm_fmadd_pd(ret, x, a0);
   return ret;
}

请注意，对于相同数量的展开项，如果针对特定的 x 范围对函数进行近似，例如使用最小二乘法，您可以获得更好的近似值。

所有这些方法都适用于非常有限的 x 范围，但连续导数在某些情况下可能很重要。

有一个技巧可以在 非常宽的范围 但具有明显的 分段线性 区域中近似指数。它基于将整数重新解释为浮点数。为了更准确的描述，我推荐这个参考：

Piecewise linear approximation to exponential and logarithm

A Fast, Compact Approximation of the Exponential Function

这种方法的可能实现方式：

__m128d fastExp3(__m128d x)
{
   const __m128d a = _mm_set1_pd(1.0 / M_LN2);
   const __m128d b = _mm_set1_pd(3 * 1024.0 - 1.05);
   __m128d t = _mm_fmadd_pd(x, a, b);
   return _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(t), 11));
}

尽管此方法简单且范围广泛 x，但在数学中使用时要小心。在小范围内，它会给出分段近似，这会破坏敏感算法，尤其是那些使用微分的算法。

要比较不同方法的准确性，请查看图形。第一张图是针对 x = [0..1) 范围制作的。如您所见，这种情况下的最佳近似值由方法 fastExp2(x) 给出，稍差但可以接受的是 fastExp1(x)。 fastExp3(x) 提供的最差近似值 - 分段结构很明显，存在一阶导数的不连续性。

在范围 x = [0..10) fastExp3(x) 方法提供了最好的近似值，更差的是 fastExp1(x) 给出的近似值 - 在相同的计算次数下，它提供了比 [=16 更多的顺序=].

下一步是提高 fastExp3(x) 算法的准确性。最简单的显着提高精度的方法是使用相等exp(x) = exp(x/2)/exp(-x/2)虽然增加了计算量，但是大大减少了除法时由于相互误差补偿造成的误差。

__m128d fastExp5(__m128d x)
{
   const __m128d ap = _mm_set1_pd(0.5 / M_LN2);
   const __m128d an = _mm_set1_pd(-0.5 / M_LN2);
   const __m128d b = _mm_set1_pd(3 * 1024.0 - 1.05);
   __m128d tp = _mm_fmadd_pd(x, ap, b);
   __m128d tn = _mm_fmadd_pd(x, an, b);
   tp = _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(tp), 11));
   tn = _mm_castsi128_pd(_mm_slli_epi64(_mm_castpd_si128(tn), 11));
   return _mm_div_pd(tp, tn);
}

通过使用等式 exp(x+dx) = exp(x) *exp(dx) 组合 fastExp1(x) 或 fastExp2(x) 和 fastExp3(x) 算法的方法，可以实现更高的准确性。如上所示，第一个乘数可以用类似fastExp3(x)的方法计算，对于第二个乘数可以使用fastExp1(x)或fastExp2(x)方法。在这种情况下找到最佳解决方案是一项艰巨的任务，我建议查看答案中提出的库中的实现。

Answer 2

有几个库提供向量化指数，或多或少的准确性。

SVML，随英特尔编译器提供（它也提供内在函数，所以如果你有许可证，你可以使用它们），具有不同级别的精度（和速度）
您提到 IPP，同样来自 Intel，它也提供一些功能
MKL 也为这个计算提供了一些接口（对于这个，修复 ISA 可以通过宏来完成，例如，如果你需要再现性或精度）
fmath 是另一种选择，您可以从矢量化 exp 中提取代码以将其集成到循环中。

根据经验，所有这些都比自定义 padde 近似更快更精确（甚至没有谈论不稳定的泰勒展开，它会很快给你负数）。

对于 SVML、IPP 和 MKL，我会检查哪一个更好：从循环内部调用还是通过对整个数组的一次调用调用 exp（因为库可以使用 AVX512 而不仅仅是 SSE2）。

Answer 3

没有 exp 的 SSE2 实现，所以如果您不想按照上面的建议自行实现，一种选择是在某些支持 ERI（指数和倒数指令）的硬件上使用 AVX512 指令。参见 https://en.wikipedia.org/wiki/AVX-512#New_instructions_in_AVX-512_exponential_and_reciprocal

我认为目前您只能使用 Xeon phi（正如 Peter Cordes 所指出的 - 我确实发现了一个关于它在 Skylake 和 Cannonlake 上的说法，但无法证实），同时请记住该代码在其他架构上根本不起作用（即会崩溃）。

你如何用 SSE2 处理 exp()？

How do you process exp() with SSE2?

c++

simd

intrinsics

sse2

exp