c - 仅使用单精度 float 逼近 [0,pi] 上的余弦

Question

我目前正在研究余弦的近似值。由于最终目标设备是使用 32 位浮点 ALU/LU 的自行开发，并且有专门的 C 编译器，因此我无法使用 C 库数学函数(cosf，...)。我的目标是编写在准确性和指令/周期数量方面不同的各种方法。
我已经尝试了很多不同的近似算法，从 fdlibm、泰勒展开、pade 近似、使用枫树的 remez 算法等开始......
但是，一旦我仅使用浮点精度来实现它们，就会显着降低精度。并且可以肯定:我知道 double ，更高的精度完全没有问题......
现在，我有一些近似值，在 pi/2(出现最大误差的范围)附近精确到几千 ulp，我觉得我受到单精度转换的限制。
为了解决主题参数减少:输入是弧度。我假设参数减少会由于除法/乘法而导致更多的精度损失......因为我的整体输入范围只有 0..pi，我决定将参数减少到 0..pi/2。
因此我的问题是:有没有人知道高精度(在最好的情况下效率高)余弦函数的单精度近似值？是否有任何算法可以优化单精度的近似值？你知道内置的 cosf 函数是否在内部计算单精度 double 的值吗？
~

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

当然可以仅使用 native 精度操作来计算 [0, π] 上的余弦，并且具有任何所需的误差界限 >= 0.5 ulp。然而，目标越接近正确舍入的函数，运行时需要的前期设计工作和计算工作就越多。
先验函数的实现通常包括参数缩减、核心近似、最终修复以抵消参数缩减。在参数减少涉及减法的情况下，需要通过显式或隐式使用更高的精度来避免灾难性取消。隐式技术可以设计为仅依赖于 native 精度计算，例如将 π 之类的常数拆分为未计算的和，例如 1.57079637e+0f - 4.37113883e-8f使用 IEEE-754 时 binary32 (单精度)。
当硬件提供融合乘加 (FMA) 运算时，通过 native 精度计算实现高精度要容易得多。 OP 没有指定他们的目标平台是否提供此操作，因此我将首先展示一种非常简单的方法，该方法仅依靠乘法和加法来提供中等精度(最大误差 < 5 ulps)。我假设硬件符合 IEEE-754 标准，并假设 float映射到 IEEE-754 binary32格式。
以下内容基于 Colin Wallace 题为“使用切比雪夫多项式逼近 sin(x) 到 5 ULP”的博客文章，在撰写本文时该文章无法在线获取。我最初检索它here并且 Google 目前保留了一个缓存副本 here .他们建议通过使用 sin(x)/(x*(x²-π²)) 的 x² 中的多项式，然后将其乘以 x*(x²-π²) 来逼近 [-π, π] 上的正弦。更准确地计算 a²-b² 的标准技巧是将其重写为 (a-b) * (a+b)。将 π 表示为两个浮点数 pi_high 和 pi_low 的未计算和避免了减法过程中的灾难性抵消，这将计算 x²-π² 变成了 ((x - pi_hi) - pi_lo) * ((x + pi_hi) + pi_lo) .
理想情况下，多项式核心近似应使用极小极大近似，即 分钟 模仿 最大 错误。我在这里已经这样做了。为此可以使用各种标准工具，例如 Maple 或 Mathematics，或者根据 Remez 算法创建自己的代码。
对于 [0, PI] 上的余弦计算，我们可以利用 cos (t) = sin (π/2 - t) 这一事实。将 x = (π/2 - t) 代入 x * (x - π/2) * (x + π/2) 产生 (π/2 - t) * (3π/2 - t) * (-π/2) - 吨)。常量可以像以前一样分为高低部分(或头和尾，使用另一种常见的习语)。

/* 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;
}

下面我展示了一种经典方法，它首先在记录象限时将参数减少到 [-π/4, π/4]。然后象限告诉我们是否需要在这个主要近似区间上计算正弦或余弦的多项式近似，以及我们是否需要翻转最终结果的符号。此代码假定目标平台支持 IEEE-754 指定的 FMA 操作，并通过标准 C 函数 fmaf() 进行映射。对于单精度。
代码很简单，除了用于计算象限的舍入模式为最近或偶数的 float-to-int 转换，这是通过“魔数(Magic Number)加法”方法执行的，并结合 2/的乘法π(相当于除以 π/2)。最大误差小于 1.5 ulps。

/* 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.451367 ulps，即保持在 1.5 ulps 以下。