linux - 为什么在 Linux/ARMv7 上编译我的代码时 sinf() 这么慢？

Question

当我在 ARMv7 上编译我的程序时，与在 Windows 上获得的相比，具有数学 sin() 函数的相同 C 代码运行速度要慢得多。我正在使用 -O2 -Wall -Wextra -mfpu=neon -mtune=cortex-a9 -march=armv7 -std=c++11 编译，我的 gcc 是 gcc (Ubuntu/Linaro 4.8.2-19ubuntu1) 4.8.2.

我不认为这只是 sin() 对于实时操作来说不是那么快，而且我知道使用查找表可以更好地折衷更快的 sin 函数，但我在这里遇到的可能是编译器中的异常或错误，因为运行 sin() 函数确实需要很长时间。

我的程序在启动时创建了几个波表，虽然它在 Windows 上几乎立即启动，但在 Linux/ARM 上启动大约需要 25-30 秒...

这里有一些代码显示了在哪里使用了 sinf() 函数来减慢一切。

for (int n = 0; n < 73; ++n)
{
    // Max number of harmonics
    int hrm = int(16000.f / twf[n]);

    // Set vectors
    basic_wf.assign(wavelength[n], 0);

    for (int i = 0; i < wavelength[n]; ++i)
    {
        // Add harmonics
        for (int h = 1; h < hrm; ++h)
        {
            const float harm = 0.14f * (sinf((float)i * FACTOR * twf[n] * (float)h) / (float)h);
            if (h % 2 == 0) basic_wf[i] -= harm;    // add even negative harmonic
            else basic_wf[i] += harm;               // add odd positive harmonic
        }
    }
}

在这里，我用锯齿波填充了 73 个表格，为每个频率添加了所需数量的谐波。音符的音高越低，泛音的数量就越高(实际的 sin() 计算)。这在 Windows 上几乎可以立即运行……在我的 Linux 机器上需要一生。

Answer 1

代码表明，您在评论中的分析证实了 sinf() 的参数量级可能会变得相当大，肯定会达到几千。三角函数的公共(public)库实现中使用的精确参数缩减可能需要大量计算，因此对于大参数来说速度很慢，尤其是当硬件平台不支持融合乘加运算时。这可能是您观察到的低 sinf() 性能的一个促成因素。

您在评论中提到 sinf() 的操作数包括一个因子 π。这表明您实际上想要使用 sinpif()，其中 sinpi(x) = sin(x*π)。 sinpi 函数在 IEEE-754 (2008) 浮点标准中引入，但尚未进入语言标准。然而，一些工具链将其作为扩展提供。 sinpi() 的优点是无论参数的大小如何，它只需要非常简单的参数约简，这可以大大减少执行时间。这会提高性能。由于乘以 π 是隐式的，因此与使用 sinf() 的离散方法相比，它还可以提供更高的精度。

我在下面展示了 sinpif() 的示例性 C99 实现。请注意，此代码在很大程度上依赖于标准数学函数 fmaf() 来实现高处理速度和出色的准确性。如果您的 CPU 不支持融合乘加 (FMA) 运算的硬件，则此函数的执行速度会非常慢，因为正确模拟 fmaf() 并非易事。由于代码是以模块化方式编写的，您可能希望将编译器配置为应用最大数量的函数内联，或者向所有组成函数添加适当的内联属性。

由于您指出您的硬件平台不提供对 FMA 的原生支持，您可以将每个 fmaf(a,b,c) 替换为 (a*b+c)，在准确性上有所损失。根据我的测试，最大 ulp 误差增加到 1.71364 ulp。这仍然非常好，但是 my_sinf() 在那种情况下不再忠实地四舍五入，这通常被认为是一个理想的属性。

/* Argument reduction for sinpi, cospi, sincospi. Reduces to [-0.25, +0.25] */
float trig_red_pi_f (float a, int *i)
{
    float r;
    r = rintf (a + a);
    *i = (int)r;
    r = a - 0.5f * r;
    return r;
}

/* Approximate cos(pi*x) for x in [-0.25,0.25]. Maximum ulp error = 0.87440 */
float cospif_poly (float s)
{
    float r;
    r =              0x1.d98dcep-3f;   //  2.31227502e-1f
    r = fmaf (r, s, -0x1.55c4e8p+0f);  // -1.33503580e+0f
    r = fmaf (r, s,  0x1.03c1d4p+2f);  //  4.05870533e+0f
    r = fmaf (r, s, -0x1.3bd3ccp+2f);  // -4.93480206e+0f
    r = fmaf (r, s,  0x1.000000p+0f);  //  1.00000000e+0f
    return r;
}

/* Approximate sin(pi*x) for x in [-0.25,0.25]. Maximum ulp error = 0.96441 */
float sinpif_poly (float a, float s)
{
    float r;
    r =             -0x1.2dc6f8p-1f;   // -5.89408636e-1f
    r = fmaf (r, s,  0x1.46602ep+1f);  //  2.54981017e+0f
    r = fmaf (r, s, -0x1.4abbc0p+2f);  // -5.16770935e+0f
    r = r * s;
    r = fmaf (r, a, -0x1.777a5cp-24f * a); // PI_lo // -8.74227766e-8f
    r = fmaf (a, 0x1.921fb6p+1f, r);       // PI_hi //  3.14159274e+0f
    return r;
}

/* Compute sin(pi*x) and cos(pi*x) based on quadrant */
float sinpif_cospif_core (float a, int i)
{
    float r, s;
    s = a * a;
    r = (i & 1) ? cospif_poly (s) : sinpif_poly (a, s);
    if (i & 2) {
        r = 0.0f - r; // don't change "sign" of NaNs or create negative zeros
    }
    return r;
}

/* maximum ulp error = 0.96411 */
float my_sinpif (float a)
{
    float r;
    int i;
    r = trig_red_pi_f (a, &i);
    r = sinpif_cospif_core (r, i);
    /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */
    r = (a == truncf (a)) ? (a * 0.0f) : r;
    return r;
}