用单精度浮点逼近[0,pi]上的余弦

Question

我目前正在研究余弦的近似值。由于最终的目标设备是一个使用32位浮点ALU/LU的自开发设备，并且有一个专门的C编译器，所以我不能使用C库的数学函数(cosf，...)。我的目标是编写在精度和指令/周期数量方面不同的各种方法。

我已经尝试了很多不同的近似算法，从fdlibm开始，taylor展开，pade近似，remezhtml" target="_blank">算法使用maple等等....

但是一旦我只使用浮点精度来实现它们，精度就会有很大的损失。而且要确定：我知道有了双倍精度，高得多的精度根本没有问题...

现在，我有一些近似，精确到几千个ulp在pi/2附近（最大误差发生的范围），我觉得我受到了单一精度转换的限制。

float ua_cos_v2(float x)
{
    float output;
    float myPi = 3.1415927410125732421875f;
    if (x < 0) x = -x;
    int quad = (int32_t)(x*0.63661977236f);//quad = x/(pi/2) = x*2/pi
    if (x<1.58f && x> 1.57f) //exclude approximation around pi/2
    {
        output = -(x - 1.57079637050628662109375f) - 2.0e-12f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 0.16666667163372039794921875f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f) + 2.0e-13f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)+ 0.000198412701138295233249664306640625f*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f)*(x - 1.57079637050628662109375f);
        output -= 4.37E-08f;
    }
    else {
        float param_x;
        int param_quad = -1;
        switch (quad)
        {
        case 0:
            param_x = x;
            break;
        case 1:
            param_x = myPi - x;
            param_quad = 1;
            break;
        case 2:
            param_x = x - myPi;
            break;
        case 3:
            param_x = 2 * myPi - x;
            break;
        }
        float c1 = 1.0f,
            c2 = -0.5f,
            c3 = 0.0416666679084300994873046875f,
            c4 = -0.001388888922519981861114501953125f,
            c5 = 0.00002480158218531869351863861083984375f,
            c6 = -2.75569362884198199026286602020263671875E-7f,
            c7 = 2.08583283978214240050874650478363037109375E-9f,
            c8 = -1.10807162057025010426514199934899806976318359375E-11f;
        float _x2 = param_x * param_x;
        output = c1 + _x2*(c2 + _x2*(c3 + _x2*(c4 + _x2*(c5 + _x2*(c6 + _x2*(c7 
        + _x2* c8))))));
        if (param_quad == 1 || param_quad == 0)
            output = -output;
    }
    return output;
}

如果我忘记了任何信息，请不要犹豫询问！

提前致谢

Answer 1

当然，只要使用本地精度运算，就可以在[0,π]上以任何期望的误差范围>=0.5ulp计算余弦。然而，目标越接近一个正确的舍入函数，就需要更多的前期设计工作和运行时的计算工作。

超越函数的实现通常包括变元约简、核心近似和最终修正来抵消变元约简。在参数减少涉及减法的情况下，需要通过显式或隐式地使用更高的精度来避免灾难性的抵消。隐式技术可以html" target="_blank">设计成只依赖于本地精度计算，例如，在使用IEEE-754binary32（单精度）时，将像π这样的常数拆分为一个未求和，如1.57079637 e+0f-4.37113883 e-8f。

当硬件提供融合乘加(FMA)操作时，用本机精度计算实现高精度要容易得多。OP没有具体说明他们的目标平台是否提供这种操作，所以我将首先展示一种非常简单的方法，它只依赖于乘法和加法，提供中等精度（最大误差<5ulps)。我假设硬件符合IEEE-754标准，并假设float映射到IEEE-754binary32格式。

以下是基于Colin Wallace的一篇博客文章，标题是“用Chebyshev多项式逼近sin(x)到5 ULP”，在撰写本文时，这篇文章还没有在网上找到。我最初是在这里检索的，Google现在在这里保留了一个缓存副本。他们建议通过使用sin(x)/(x*(x2--π2))中的一个多项式来逼近[-π,π]上的正弦，然后将其乘以x*(x2--π2))。要更准确地计算Aâ-Bâ，一个标准技巧是将其改写为(a-b)*(A+b)。将π表示为两个浮点数pi_high和pi_low的未求值之和，可以避免减法过程中的灾难性抵消，从而将计算xè-πé变为((x-pi_hi)-pi_lo)*((x+pi_hi)+pi_lo)。

多项式核近似最好使用极小极大近似，使最大误差最小。我在这里已经这么做了。可以使用各种标准工具，如Maple或数学，或者根据Remez算法创建自己的代码。

对于[0,PI]上的余弦计算，我们可以利用cos(t)=sin(π/2-t)这一事实。将x=(π/2-t)代入x*(x-π/2)*(x+π/2)得到(π/2-t)*(3π/2-t)*(-π/2-t)。常量可以像前面一样分为高和低两部分（或者用另一个常见的习惯用法来说，头和尾）。

/* Approximate cosine on [0, PI] with maximum error of 4.704174 ulp */
float cosine (float x)
{
    const float half_pi_hi       =  1.57079637e+0f; //  0x1.921fb6p+0
    const float half_pi_lo       = -4.37113883e-8f; // -0x1.777a5cp-25
    const float three_half_pi_hi =  4.71238899e+0f; //  0x1.2d97c8p+2
    const float three_half_pi_lo = -1.19248806e-8f; // -0x1.99bc5cp-27
    float p, s, hpmx, thpmx, nhpmx;

    /* cos(x) = sin (pi/2 - x) = sin (hpmx) */
    hpmx = (half_pi_hi - x) + half_pi_lo;               // pi/2-x
    thpmx = (three_half_pi_hi - x) + three_half_pi_lo;  // 3*pi/2 - x
    nhpmx = (-half_pi_hi - x) - half_pi_lo;             // -pi/2 - x

    /* P(hpmx*hpmx) ~= sin (hpmx) / (hpmx * (hpmx * hpmx - pi * pi)) */
    s = hpmx * hpmx;
    p =         1.32729383e-10f;
    p = p * s - 2.33177868e-8f;
    p = p * s + 2.52223435e-6f;
    p = p * s - 1.73503853e-4f;
    p = p * s + 6.62087463e-3f;
    p = p * s - 1.01321176e-1f;
    return hpmx * nhpmx * thpmx * p;
}

代码简单明了，除了用于计算象限的从浮点到int的转换，该转换采用舍入模式到最近或偶，该转换通过“幻数加法”方法执行，并与2/π的乘法相结合（等效于π/2)。最大误差小于1.5ulps。

/* compute cosine on [0, PI] with maximum error of 1.429027 ulp */
float my_cosf (float a)
{
    const float half_pi_hi =  1.57079637e+0f; //  0x1.921fb6p+0
    const float half_pi_lo = -4.37113883e-8f; // -0x1.777a5cp-25
    float c, j, r, s, sa, t;
    int i;

    /* subtract closest multiple of pi/2 giving reduced argument and quadrant */
    j = fmaf (a, 6.36619747e-1f, 12582912.f) - 12582912.f; // 2/pi, 1.5 * 2**23
    a = fmaf (j, -half_pi_hi, a);
    a = fmaf (j, -half_pi_lo, a);

    /* phase shift of pi/2 (one quadrant) for cosine */
    i = (int)j;
    i = i + 1;

    sa = a * a;
    /* Approximate cosine on [-PI/4,+PI/4] with maximum error of 0.87444 ulp */
    c =               2.44677067e-5f;  //  0x1.9a8000p-16
    c = fmaf (c, sa, -1.38877297e-3f); // -0x1.6c0efap-10
    c = fmaf (c, sa,  4.16666567e-2f); //  0x1.555550p-5
    c = fmaf (c, sa, -5.00000000e-1f); // -0x1.000000p-1
    c = fmaf (c, sa,  1.00000000e+0f); //  1.00000000p+0
    /* Approximate sine on [-PI/4,+PI/4] with maximum error of 0.64196 ulp */
    s =               2.86567956e-6f;  //  0x1.80a000p-19
    s = fmaf (s, sa, -1.98559923e-4f); // -0x1.a0690cp-13
    s = fmaf (s, sa,  8.33338592e-3f); //  0x1.111182p-7
    s = fmaf (s, sa, -1.66666672e-1f); // -0x1.555556p-3
    t = a * sa;
    s = fmaf (s, t, a);

    /* select sine approximation or cosine approximation based on quadrant */
    r = (i & 1) ? c : s;
    /* adjust sign based on quadrant */
    r = (i & 2) ? (0.0f - r) : r;

    return r;
}

事实证明，在这种特殊情况下，FMA的使用在准确性方面只提供了一个微小的好处。如果将对FMAF(a，b，c)的调用替换为((a)*(b)+(c))，则最大错误最小地增加到1.451367ulps，即保持在1.5ulps以下。