c - 为等间距角度快速准确地迭代生成正弦和余弦

Question

在某些应用中，需要多个角度的正弦和余弦，其中角度是通过将相等大小的增量 incr 重复添加到起始值 base 得出的。出于性能原因，而不是调用 sin()，cos() 标准数学库函数(或者可能是非标准的 sincos()函数)对于每个生成的角度，只计算一次 sin(base) 和 cos(base)，然后通过应用导出所有其他正弦和余弦是非常有利的angle-sum formulas :

sin(base+incr) = cos(incr) · sin(base) + sin(incr) · cos(base)
cos(base+incr) = cos(incr) · cos(base) - sin(incr) · sin(base)

这只需要对比例因子 sin(incr) 和 cos(incr) 进行一次性预计算，而不管执行了多少次迭代。

这种方法有几个问题。如果增量很小，cos(incr) 将是一个接近于 1 的数字，当以有限精度浮点格式计算时，通过隐式减法抵消导致精度损失。此外，由于计算未以有利的数值形式 sin(base+incr) = sin(base) + adjust 进行计算，因此会产生比必要更多的舍入误差，其中计算量 adjust 的量级明显小于 sin(base)(类似于余弦)。

由于一个人通常应用数十到数百次迭代步骤，因此这些错误会累积。如何以最有利于保持高精度的方式构建迭代计算？如果 fused multiply-add 应该对算法进行哪些更改操作 (FMA) 可用，它通过标准数学函数 fma() 和 fmaf()?

公开

Answer 1

正弦的 half-angle formula 应用允许解决问题中提到的两个影响准确性的问题:

sin(incr/2) = √((1-cos(incr))/2) ⇒
sin²(incr/2) = (1-cos(incr))/2 ⇔
2·sin²(incr/2) = 1-cos(incr) ⇔
1-2·sin²(incr/2) = cos(incr)

将其代入原始公式会得到此中间表示:

sin(base+incr) = (1 - 2·sin²(incr/2)) · sin(base) + sin(incr) · cos(base)
cos(base+incr) = (1 - 2·sin²(incr/2)) · cos(base) - sin(incr) · sin(base)

通过对项进行简单的重新排序，可以得到所需的公式形式:

sin(base+incr) = sin(base) + (sin(incr) · cos(base) - 2·sin²(incr/2) · sin(base))
cos(base+incr) = cos(base) - (2·sin²(incr/2) · cos(base) + sin(incr) · sin(base))

和原来的公式一样，这只需要一次性预计算两个比例因子，即2·sin²(incr/2)和sin(incr) .对于小增量，两者都很小:保留完整的精度。

关于如何将 FMA 应用于此计算，有两种选择。可以通过废除使用一次调整的方法来最小化操作次数，而是使用两次，希望 FMA 操作的减少舍入误差(一次舍入用于两次操作)将补偿精度损失:

sin(base+incr) = fma (-2·sin²(incr/2), sin(base), fma ( sin(incr), cos(base), sin(base)))
cos(base+incr) = fma (-2·sin²(incr/2), cos(base), fma (-sin(incr), sin(base), cos(base)))

另一种选择是将单个 FMA 应用于改进后的公式，尽管尚不清楚应将两个乘法中的哪一个映射到 FMA 内的未舍入乘法:

sin(base+incr) = sin(base) + fma (sin(incr), cos(base), -2·sin²(incr/2) · sin(base))
cos(base+incr) = cos(base) - fma (sin(incr), sin(base), 2·sin²(incr/2) · cos(base))

下面的脚手架通过生成许多 (base, incr) 对来评估上面讨论的每个计算备选方案，然后对它们中的每一个迭代一定数量的步骤同时收集生成的所有正弦和余弦值的误差。由此它为每个测试用例计算一个 root-mean square error，分别为正弦、余弦。最后报告在所有生成的测试用例中观察到的最大 RMS 误差。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NAIVE    (1)
#define ROBUST   (2)
#define FAST     (3)
#define ACCURATE (4)
#define MODE (ACCURATE)

// Fixes via: Greg Rose, KISS: A Bit Too Simple. http://eprint.iacr.org/2011/007
static unsigned int z=362436069,w=521288629,jsr=362436069,jcong=123456789;
#define znew (z=36969*(z&0xffff)+(z>>16))
#define wnew (w=18000*(w&0xffff)+(w>>16))
#define MWC  ((znew<<16)+wnew)
#define SHR3 (jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5)) /* 2^32-1 */
#define CONG (jcong=69069*jcong+13579)                     /* 2^32 */
#define KISS ((MWC^CONG)+SHR3)

int main (void)
{
    double sumerrsqS, sumerrsqC, rmsS, rmsC, maxrmsS = 0, maxrmsC = 0;
    double refS, refC, errS, errC;
    float base, incr, s0, c0, s1, c1, tt;
    int count, i;

    const int N = 100;  // # rotation steps per test case
    count = 2000000;    // # test cases (a pair of base and increment values)

#if MODE == NAIVE
    printf ("testing: NAIVE (without FMA)\n");
#elif MODE == FAST
    printf ("testing: FAST (without FMA)\n");
#elif MODE == ACCURATE
    printf ("testing: ACCURATE (with FMA)\n");
#elif MODE == ROBUST
    printf ("testing: ROBUST (with FMA)\n");
#else
#error unsupported MODE
#endif // MODE

    do {
        /* generate test case */
        base = (float)(KISS * 1.21e-10);      // starting angle, < 30 degrees
        incr = (float)(KISS * 2.43e-10 / N);  // increment, < 60/n degrees

        /* set up rotation parameters */
        s1 = sinf (incr);
#if MODE == NAIVE
        c1 = cosf (incr);
#else
        tt = sinf (incr * 0.5f);
        c1 = 2.0f * tt * tt;
#endif // MODE
        sumerrsqS = 0;
        sumerrsqC = 0;

        s0 = sinf (base); // initial sine
        c0 = cosf (base); // initial cosine

        /* run test case through N rotation steps */
        i = 0;
        do {         

            tt = s0; // old sine
#if MODE == NAIVE
            /* least accurate, 6 FP ops */
            s0 = c1 * tt + s1 * c0; // new sine
            c0 = c1 * c0 - s1 * tt; // new cosine
#elif MODE == ROBUST
            /* very accurate, 8 FP ops */
            s0 = ( s1 * c0 - c1 * tt) + tt; // new sine
            c0 = (-s1 * tt - c1 * c0) + c0; // new cosine
#elif MODE == FAST
            /* accurate and fast, 4 FP ops */
            s0 = fmaf (-c1, tt, fmaf ( s1, c0, tt)); // new sine
            c0 = fmaf (-c1, c0, fmaf (-s1, tt, c0)); // new cosine
#elif MODE == ACCURATE
            /* most accurate, 6 FP ops */
            s0 = tt + fmaf (s1, c0, -c1 * tt); // new sine
            c0 = c0 - fmaf (s1, tt,  c1 * c0); // new cosine
#endif // MODE
            i++;

            refS = sin (fma ((double)i, (double)incr, (double)base));
            refC = cos (fma ((double)i, (double)incr, (double)base));
            errS = ((double)s0 - refS) / refS;
            errC = ((double)c0 - refC) / refC;
            sumerrsqS = fma (errS, errS, sumerrsqS);
            sumerrsqC = fma (errC, errC, sumerrsqC);
        } while (i < N);

        rmsS = sqrt (sumerrsqS / N);
        rmsC = sqrt (sumerrsqC / N);
        if (rmsS > maxrmsS) maxrmsS = rmsS;
        if (rmsC > maxrmsC) maxrmsC = rmsC;

    } while (--count);

    printf ("max rms error sin = % 16.9e\n", maxrmsS);
    printf ("max rms error cos = % 16.9e\n", maxrmsC);

    return EXIT_SUCCESS;
}

测试支架的输出表明，最快的基于 FMA 的替代方案优于问题中的朴素方法，而更准确的基于 FMA 的替代方案是所考虑替代方案中最准确的:

testing: NAIVE (without FMA)
max rms error sin =  4.837386842e-006
max rms error cos =  6.884047862e-006

testing: ROBUST (without FMA)
max rms error sin =  3.330292645e-006
max rms error cos =  4.297631502e-006

testing: FAST (with FMA)
max rms error sin =  3.532624939e-006
max rms error cos =  4.763623188e-006

testing: ACCURATE (with FMA)
max rms error sin =  3.330292645e-006
max rms error cos =  4.104813533e-006