如何使用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);
}

我会说它运作得很好.但是,找不到expSSE2的任何功能,来完成操作链.

读这个,似乎我需要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 范围内工作，但具有在某些情况下可能很重要的连续导数。

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

指数和对数的分段线性逼近

指数函数的快速、紧凑逼近

这种方法的可能实现：

__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);
}

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

Answer 2

有几个库提供矢量化指数,具有或多或少的准确性.

• SVML,随英特尔编译器提供(它也提供内在函数,所以如果你有许可证,你可以使用它们),具有不同的精度(和速度)
• 你提到了同样来自英特尔的IPP,它也提供了一些功能
• MKL还为此计算提供了一些接口(对于这一点,修复ISA可以通过宏完成,例如,如果您需要可重复性或精度)
• fmath是另一种选择,你可以从vectorized exp中删除代码,将它集成到你的循环中.

根据经验,所有这些都比自定义padde近似更快更精确(甚至没有谈论不稳定的泰勒扩展,这会非常快地给你负数).

对于SVML,IPP和MKL,我会检查哪一个更好:从你的循环内部调用或调用exp来调用你的完整数组(因为库可以使用AVX512而不仅仅是SSE2).