仅使用单精度浮点在 [0,pi] 上逼近余弦答案

【问题标题】：Approximating cosine on [0,pi] using only single precision floating point仅使用单精度浮点在 [0,pi] 上逼近余弦
【发布时间】：2021-01-03 05:11:00
【问题描述】：

我目前正在研究余弦的近似值。由于最终目标设备是使用 32 位浮点 ALU / LU 的自行开发，并且有专门的 C 编译器，因此我无法使用 c 库数学函数（cosf，...）。我的目标是编写在准确性和指令/周期数方面不同的各种方法。

我已经尝试了很多不同的逼近算法，从 fdlibm、taylor 展开、pade 逼近、remez 算法使用 maple 等等......

但是，一旦我只使用浮点精度来实现它们，精度就会大大降低。并且可以肯定：我知道使用双精度，更高的精度完全没有问题......

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

为了解决主题参数减少：输入以弧度为单位。我假设参数减少会由于除法/乘法而导致更多的精度损失......因为我的整体输入范围只有 0..pi，我决定将参数减少到 0..pi/2。

因此我的问题是：有没有人知道高精度的余弦函数的单精度近似（最好的情况是高效率）？是否有任何算法可以优化单精度近似值？你知道内置的 cosf 函数是否在内部以单精度或双精度计算值？ ~

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

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

提前致谢

【问题讨论】：

您需要什么精度？请显示精度不足的近似代码。也许有一些方法可以提高精度。（不看代码我们无法判断。）请edit你的问题添加此信息，不要使用cmets来回答。
您的实际输入是以弧度为单位的，还是您真的想为0<=x<=1 计算cos(x*pi)？无论如何，在应用任何类型的多项式近似之前，您应该将输入范围缩小到[-pi/4, pi/4] 并使用诸如cos(x+pi/2) = -sin(x) 之类的标识。
对于接近 π/2 的 x，也就是你所说的最大误差，cos(x) 在 π/2−x 附近。这意味着用多项式近似它很容易。具体来说，您应该使用 y=π/2−x，然后这种情况下的特定多项式是 y，但即使是某种形式的更一般的多项式，例如 y+c3•y^3+c5•y^5+… 也会有由于高阶项实际上为零，因此计算误差很小。发生错误的地方是计算 y=π/2−x。如果以高精度完成此操作，则 y 的结果在 π/2 附近的 ULP 的一小部分内是准确的。如果使用float 精度完成，则错误很大。
对于这种特定情况，您可以考虑将 π/2 分为两部分存储。第一个 p0 是 π/2 舍入到 float。第二个，p1，是 π/2−p0（预先计算，结果写入源代码）。然后 π/2−x 可以在float 精度和p0-x+p1 中精确计算。当 x 是最接近 π/2 的float 时，这会产生大约 ⅓ ULP 的误差。除此之外，我们还需要查看您正在使用的代码。
“高精度（最好的情况下是高效率）”——>哪个更重要？我建议最快，只要准确度为

标签： c floating-point trigonometry approximation single-precision

【解决方案1】：

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

超越函数实现通常包括参数缩减、核心逼近、最终修复以抵消参数缩减。在参数减少涉及减法的情况下，需要通过显式或隐式使用更高的精度来避免灾难性的取消。隐式技术可以设计为仅依赖本机精度计算，例如，在使用 IEEE-754 binary32（单精度）时，将像 π 这样的常数拆分为未计算的和，如 1.57079637e+0f - 4.37113883e-8f。

当硬件提供融合乘加 (FMA) 操作时，使用本机精度计算实现高精度会容易得多。 OP 没有说明他们的目标平台是否提供此操作，所以我将首先展示一个非常简单的方法，仅依靠乘法和加法来提供中等精度（最大误差 float 映射到 IEEE-754 binary32 格式。

以下内容基于 Colin Wallace 题为“使用 Chebyshev 多项式将 sin(x) 近似为 5 ULP”的博客文章，该文章在撰写本文时尚未在线提供。我最初检索了它here，而谷歌目前保留了一个缓存副本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)。

多项式核心近似应该理想地使用一个极小极大近似，它min将max最大误差最小化。我在这里这样做了。为此可以使用 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;
}

下面我展示了一种经典方法，它首先在记录象限时将参数简化为 [-π/4, π/4]。然后象限告诉我们是否需要在这个主要近似区间上计算正弦或余弦的多项式近似，以及我们是否需要翻转最终结果的符号。此代码假定目标平台支持 IEEE-754 指定的 FMA 操作，并通过标准 C 函数 fmaf() 映射为单精度。

代码很简单，除了用于计算象限的舍入模式到最近或偶数的浮点到整数转换，这是通过“幻数加法”方法执行并与乘法相结合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 在准确性方面只提供了很小的好处。如果我用((a)*(b)+(c)) 替换对fmaf(a,b,c) 的调用，最大误差最小地增加到1.451367 ulps，也就是说，它保持在1.5 ulps 以下。

【讨论】：

非常感谢！我实际上必须检查硬件是否提供 FMA。您的实现使用了 ericpostpischil 提到的相同想法。实际上，这些几行代码如何产生这样的结果给我留下了深刻的印象。非常有趣，但又如此简单。我已经测试了您的第一种方法，它适用于大约 30 条指令，非常棒。我现在将实施您的第二个解决方案。 fma 的微小好处很合适，因为我没有可用的标准数学库.... 每个值 (a)、(b) 周围的括号是否有特定原因
@DexterS 重新括号：只需从宏定义中剪切和粘贴：#define fmaf(a,b,c) ((a)*(b)+(c))。自 1970 年代以来，用于提高中间计算有效精度的未求值和和拆分常数技术就已经存在，并且早于我的专业工作开始。原作者是 Kahan、Dekker、Cody/Waite。
刚刚测试了你的第二种方法，它工作得很好。我会接受你的回答，谢谢！如果我不想使用宏，我可以将其硬编码为例如j = ((a)* (6.36619747e-1f) +(12582912.f)) - 12582912.f; ，对吧？
@DexterS 绝对。我不建议使用宏来覆盖生产代码中的标准 C 函数。由于您是新提问者，我会提到，在接受答案之前至少等待 24 小时被认为是一种很好的形式，以便为所有时区的答案提供者提供平等的贡献机会。
谢谢。好的，我会等待接受答案。我还有两个问题，只是为了让我理解算法：您提到可以使用例如 maple minimax 轻松计算 remez 算法。所以这个想法总是首先得到双精度算法，然后在第二步通过拆分常数将方程转换为单精度？ ((a)*(b)+(c)) /= (a*b+c) 怎么样？

【解决方案2】：

我看到@njuffa 有一个很好的方法，但想提出另一种方法：

角度可能最初以度为单位，而不是弧度，并利用它。
不依赖于 float 是 IEEE。
fma 可能是weak，所以不要使用它。

使用整数数学进行范围缩减，然后通过自调整泰勒级数找到答案。

#include <assert.h>

static float my_sinf_helper(float xx, float term, unsigned n) {
  if (term + 1.0f == 1.0f) {
    return term;
  }
  return term - my_sinf_helper(xx, xx * term / ((n + 1) * (n + 2)), n + 2);
}

static float my_cosf_helper(float xx, float term, unsigned n) {
  if (term + 1.0f == 1.0f) {
    return term;
  }
  return term - xx * my_cosf_helper(xx, term / ((n + 1) * (n + 2)), n + 2);
}

// valid for [-pi/4 + pi/4]
static float my_sinf_primary(float x) {
  return x * my_sinf_helper(x * x, 1.0, 1);
}

// valid for [-pi/4 + pi/4]
static float my_cosf_primary(float x) {
  return my_cosf_helper(x * x, 1.0, 0);
}

#define MY_PIf 3.1415926535897932384626433832795f
#define D2Rf(d) ((d)*(MY_PIf/180))

float my_cosdf(float x) {
  if (x < 0) {x = -x;}
  unsigned long long ux = (unsigned long long) x;
  x -= (float) ux;
  unsigned ux_primary = ux % 360u;
  int uxq = ux_primary%90;
  if (uxq >= 45) uxq -= 90;
  x += uxq;
  switch (ux_primary/45) {
    case 7: //
    case 0: return my_cosf_primary(D2Rf(x));
    case 1: //
    case 2: return -my_sinf_primary(D2Rf(x));
    case 3: //
    case 4: return -my_cosf_primary(D2Rf(x));
    case 5: //
    case 6: return my_sinf_primary(D2Rf(x));
  }
  assert(0);
  return 0;
}

测试代码

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define DBL_FMT "%+24.17e"

typedef struct {
  double x, y0, y1, adiff;
  unsigned n;
} test;

test worst = {0};

int my_cosd_test(float x) {
  test t;
  t.x = x;
  t.y0 = cos(x*acos(-1)/180);
  t.y1 = my_cosdf(x);
  t.adiff = fabs(t.y1 - t.y0);
  if (t.adiff > worst.adiff) {
    t.n = worst.n + 1;
    printf("n:%3u x:" DBL_FMT " y0:" DBL_FMT " y1:" DBL_FMT " d:" DBL_FMT "\n", //
        t.n, t.x, t.y0, t.y1, t.adiff);
    fflush(stdout);
    worst = t;
    if (t.n > 100)
      exit(-1);
  }
  return t.adiff != 0.0;
}

float rand_float_finite(void) {
  union {
    float f;
    unsigned char uc[sizeof(float)];
  } u;
  do {
    for (size_t i = 0; i < sizeof u.uc / sizeof u.uc[0]; i++) {
      u.uc[i] = (unsigned char) rand();
    }
  } while (!isfinite(u.f) || fabs(u.f) > 5000);
  return u.f;
}

int my_cosd_tests(unsigned n) {
  my_cosd_test(0.0);
  for (unsigned i = 0; i < n; i++) {
    my_cosd_test(rand_float_finite());
  }
  return 0;
}

int main(void) {
  my_cosd_tests(1000000);
}

最坏的转换错误：+8.2e-08。最大递归深度注：6。

n: 14 x:+3.64442993164062500e+03 y0:+7.14107074054115110e-01 y1:+7.14107155799865723e-01 d:+8.17457506130381262e-08

稍后我会回顾更多。我确实看到更广泛的测试达到了大约 9e-08 最坏情况错误和x > about 1e10 的一些待定问题。

【讨论】：