您如何使用 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) 给出的近似值稍差——在相同的计算数量下，它提供的顺序比 fastExp2(x) 更多。

下一步是提高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，由 Intel 编译器提供（它也提供内部函数，所以如果你有许可证，就可以使用它们），具有不同级别的精度（和速度）
• 您提到了IPP，同样来自英特尔，它也提供了一些功能
• 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 上的声明，但无法证实），并且请记住代码在其他架构上根本不起作用（即会崩溃）。