Приближение косинуса к [0,pi] с использованием только одинарной точности с плавающей запятой

Безусловно, можно вычислить косинус на [0, π] с любой желаемой границей ошибки ›= 0,5 ulp, используя только собственные операции точности. Однако чем ближе цель к правильно округленной функции, тем больше требуется предварительных проектных и вычислительных работ во время выполнения.

Реализации трансцендентных функций обычно состоят из сокращения аргументов, основного приближения, окончательного исправления для противодействия уменьшению аргументов. В тех случаях, когда сокращение аргумента включает вычитание, необходимо избегать катастрофической отмены, явно или неявно используя более высокую точность. Неявные методы могут быть разработаны так, чтобы полагаться только на собственное вычисление точности, например, путем разделения константы, такой как π, на невычисленную сумму, такую ​​​​как 1.57079637e+0f - 4.37113883e-8f, при использовании IEEE-754 binary32 (одинарная точность).

Достичь высокой точности с помощью встроенных вычислений точности намного проще, если аппаратно реализована операция плавного умножения-сложения (FMA). OP не указал, обеспечивает ли их целевая платформа эту операцию, поэтому я сначала покажу очень простой подход, предлагающий умеренную точность (максимальная ошибка ‹ 5 ulps), основанный только на умножениях и добавлениях. Я предполагаю, что оборудование соответствует стандарту IEEE-754, и предполагаю, что float отображается в формат IEEE-754 binary32.

Нижеследующее основано на сообщении в блоге Колина Уоллеса под названием «Приближение sin(x) к 5 ULP с полиномами Чебышева», которое на момент написания недоступно в Интернете. Первоначально я получил его здесь, и в настоящее время Google сохраняет кешированную копию здесь. Они предлагают аппроксимировать синус на [-π, π] с помощью полинома от x² от sin(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] при записи квадранта. Затем квадрант сообщает нам, нужно ли нам вычислять полиномиальное приближение к синусу или косинусу на этом интервале первичного приближения, и нужно ли нам менять знак конечного результата. Этот код предполагает, что целевая платформа поддерживает операцию FMA, указанную в IEEE-754, и что она отображается через стандартную функцию 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 дает лишь незначительное преимущество с точки зрения точности. Если я заменю вызовы fmaf(a,b,c) на ((a)*(b)+(c)), максимальная ошибка увеличится минимально до 1,451367 ulps, то есть останется ниже 1,5 ulps.

Я вижу, что у @njuffa есть хороший подход, но я хочу предложить другой подход:

• Угол, скорее всего, изначально указан в градусах, а не в радианах, и воспользуйтесь этим.
• Не зависит от того, является ли float IEEE.
• fma может быть слабым и поэтому не использовать его.

Выполните сокращение диапазона, используя целочисленную математику, затем найдите ответ с помощью самонастраивающегося ряда Тейлора.

#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 и некоторых проблем TBD с x > about 1e10.