Быстрый алгоритм Arc Cos?

У меня есть моя, очень быстрая функция cos:

float sine(float x)
{
    const float B = 4/pi;
    const float C = -4/(pi*pi);

    float y = B * x + C * x * abs(x);

    //  const float Q = 0.775;
    const float P = 0.225;

    y = P * (y * abs(y) - y) + y;   // Q * y + P * y * abs(y)


    return y;
}

float cosine(float x)
{
    return sine(x + (pi / 2));
}

Но теперь, когда я просматриваю профиль, я вижу, что acos() убивает процессор. Мне не нужна интенсивная точность. Что такое быстрый способ вычисления acos (x) Спасибо.

Ответ 1

Простая кубическая аппроксимация, полином Лагранжа для x ∈ {-1, -½, 0, ½, 1} является:

double acos(x) {
   return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966;
}

Он имеет максимальную погрешность около 0,18 рад.

Ответ 2

Есть запасная память? Таблица поиска (с интерполяцией, если требуется) будет самой быстрой.

Ответ 3

nVidia имеет несколько отличных ресурсов, которые показывают, как аппроксимировать в противном случае очень дорогие математические функции, такие как: acos asin atan2 и т.д. и т.д.

Эти алгоритмы дают хорошие результаты, когда скорость выполнения важнее (в пределах разумного), чем точность. Здесь их acos функция:

// Absolute error <= 6.7e-5
float acos(float x) {
  float negate = float(x < 0);
  x = abs(x);
  float ret = -0.0187293;
  ret = ret * x;
  ret = ret + 0.0742610;
  ret = ret * x;
  ret = ret - 0.2121144;
  ret = ret * x;
  ret = ret + 1.5707288;
  ret = ret * sqrt(1.0-x);
  ret = ret - 2 * negate * ret;
  return negate * 3.14159265358979 + ret;
}

И вот результаты при вычислении acos (0.5):

nVidia:   result: 1.0471513828611643
math.h:   result: 1.0471975511965976

Это довольно близко! В зависимости от вашей требуемой степени точности это может быть хорошим вариантом для вас.

Ответ 4

У меня есть своя. Это довольно точно и быстро. Это работает с одной теоремой, построенной вокруг квадрической сходимости. Это действительно интересно, и вы можете видеть уравнение и как быстро это может привести к схождению моего естественного логарифмического приближения: https://www.desmos.com/calculator/yb04qt8jx4

Здесь мой код arccos:

function acos(x)
    local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
    local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
    local c=0.88-0.77*x c=(c+(2-a)/c)/2
    return (8*(c+(2-a)/c)-(b+(2-2*x)/b))/6
end

Многие из них являются просто квадратным корневым приближением. Он отлично работает, если вы слишком близко не используете квадратный корень из 0. Он имеет среднюю ошибку (исключая x = 0,99 до 1) 0,0003. Проблема, однако, в том, что в 0.99 она начинает гоняться, а при х = 1 разница в точности становится равной 0,05. Конечно, это можно было бы решить, выполняя больше итераций по квадратным корням (lol nope) или, только немного, например, если x > 0.99, то используйте другой набор линеаризации с квадратным корнем, но это делает код длинным и уродливым.

Если вы так не заботитесь о точности, вы можете просто сделать одну итерацию на квадратный корень, которая все равно должна держать вас где-то в диапазоне 0,0162 или что-то до точности:

function acos(x)
    local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
    local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
    local c=0.88-0.77*x c=(c+(2-a)/c)/2
    return 8/3*c-b/3
end

Если с ним все в порядке, вы можете использовать уже существующий квадратный корневой код. Он избавится от уравнения, которое немного сходит с ума при x = 1:

function acos(x)
    local a = math.sqrt(2+2*x)
    local b = math.sqrt(2-2*x)
    local c = math.sqrt(2-a)
    return 8/3*d-b/3
end

Честно говоря, если вы действительно нажали на время, помните, что вы можете линеаризовать arccos в 3.14159-1.57079x и просто делать:

function acos(x)
    return 3.14159-1.57079*x
end

В любом случае, если вы хотите увидеть список моих аппроксимационных уравнений arccos, вы можете перейти к https://www.desmos.com/calculator/tcaty2sv8l. Я знаю, что мои приближения не самые лучшие для некоторые вещи, но если вы делаете что-то, где мои аппроксимации будут полезны, используйте их, но постарайтесь дать мне кредит.

Ответ 5

Другой подход, который вы можете предпринять, - использовать сложные числа. Из de Moivre formula,

ⅈ ^x= cos (& pi;/2 * x) + ⅈ * sin (& pi;/2 * x)

Пусть & theta; = & pi;/2 * x. Тогда x = 2 & theta;/& pi, поэтому

sin (& theta;) = ℑ (ⅈ ^ ^{2 & thet;/& pi;})
cos (& theta;) = ℜ (ⅈ ^ ^{2 & theta;/& pi;})

Как вы можете вычислить степени ⅈ без sin и cos? Начните с предварительно вычисленной таблицы для степеней 2:

ⅈ ⁴= 1
ⅈ ²= -1
ⅈ ¹= ⅈ
ⅈ ^1/2= 0.7071067811865476 + 0.7071067811865475 * ⅈ
ⅈ ^1/4= 0,9238795325112867 + 0,3826834323650898 * ⅈ
ⅈ ^1/8= 0,9807852804032304 + 0,19509032201612825 * ⅈ
ⅈ ^1/16= 0.9951847266721969 + 0.0980171403295606 * ⅈ
ⅈ ^1/32= 0.9987954562051724 + 0.049067674327418015 * ⅈ
ⅈ ^1/64= 0.9996988186962042 + 0.024541228522912288 * ⅈ
ⅈ ^1/128= 0.9999247018391445 + 0.012271538285719925 * ⅈ
ⅈ ^1/256= 0.9999811752826011 + 0.006135884649154475 * ⅈ

Чтобы вычислить произвольные значения ⅈ ^x аппроксимируйте показатель как двоичную дробь и затем умножьте вместе соответствующие значения из таблицы.

Например, найти sin и cos 72 & deg; = 0,8π/2:

ⅈ ^0.8 & Около; ⅈ ^205/256 = ⅈ ^0b11001101 = ⅈ ^1/2 * ⅈ ^1/4 * ⅈ ^1/32 * ⅈ ^1/64 * ⅈ ^1/256
= 0.3078496400415349 + 0.9514350209690084 * ⅈ

sin (72 & deg;) & approx; 0,9514350209690084 ( "точное" значение составляет 0,9510565162951535).
cos (72 & deg;) & approx; 0.3078496400415349 ( "точное" значение - 0,30901699437494745).

Чтобы найти asin и acos, вы можете использовать эту таблицу с помощью метода Bisection:

Например, чтобы найти asin (0.6) (наименьший угол в 3-4-5 треугольниках):

ⅈ ⁰= 1 + 0 * ⅈ. Грех слишком мал, поэтому увеличивайте x на 1/2.
ⅈ ^1/2= 0.7071067811865476 + 0.7071067811865475 * ⅈ. Грех слишком велик, поэтому уменьшите x на 1/4.
ⅈ ^1/4= 0,9238795325112867 + 0,3826834323650898 * ⅈ. Грех слишком мал, поэтому увеличивайте x на 1/8.
ⅈ ^3/8= 0.8314696123025452 + 0,5555702330196022 * ⅈ. Грех все еще слишком мал, поэтому увеличивайте x на 1/16.
ⅈ ^7/16= 0.773010453362737 + 0,6343932841636455 * ⅈ. Грех слишком велик, поэтому уменьшите x на 1/32.
ⅈ ^13/32= 0.8032075314806449 + 0,5956993044924334 * ⅈ.

Каждый раз, когда вы увеличиваете x, умножьте на соответствующую мощность ⅈ. Каждый раз, когда вы уменьшаете x, делите на соответствующую мощность ⅈ.

Если мы остановимся здесь, получим acos (0.6) & approx; 13/32 * & pi;/2 = 0,6381360077604268 ( "Точное" значение 0,6435011087932844.)

Точность, конечно, зависит от количества итераций. Для быстрого и грязного приближения используйте 10 итераций. Для "интенсивной точности" используйте 50-60 итераций.

Ответ 6

Вы можете аппроксимировать обратный косинус с помощью полинома как предложено dan04, но многочлен является довольно плохим приближением вблизи -1 и 1, где производная от обратного косинус переходит в бесконечность. Когда вы увеличиваете степень полинома, вы быстро достигаете уменьшающейся отдачи, и по-прежнему трудно получить хорошее приближение вокруг конечных точек. В этом случае рациональная функция (фактор двух многочленов) может дать гораздо лучшее приближение.

acos(x) ≈ π/2 + (ax + bx³) / (1 + cx² + dx⁴)

где

a = -0.939115566365855
b =  0.9217841528914573
c = -1.2845906244690837
d =  0.295624144969963174

имеет максимальную абсолютную погрешность 0,017 радианов (0,96 градуса) на интервале (-1,1). Здесь график (обратный косинус в черном, кубическое полиномиальное приближение в красном, указанная выше функция в синем) для сравнения:

Приведенные выше коэффициенты выбраны так, чтобы минимизировать максимальную абсолютную погрешность по всему домену. Если вы готовы разрешить большую ошибку в конечных точках, ошибку на интервале (-0.98, 0.98) можно сделать намного меньше. Числитель степени 5 и знаменатель степени 2 примерно так же быстро, как и выше, но немного менее точный. За счет производительности вы можете повысить точность, используя полиномы более высокой степени.

Замечание о производительности: вычисление двух многочленов по-прежнему очень дешево, и вы можете использовать плавные инструкции с несколькими добавками. Деление не так уж плохо, потому что вы можете использовать аппаратное взаимное приближение и умножить. Ошибка в обратном приближении пренебрежимо мала по сравнению с ошибкой в приближении acos. На 2,6 ГГц Skylake i7 это приближение может выполнять около 8 обратных косинусов каждые 6 циклов с использованием AVX. (То есть пропускная способность, латентность длиннее 6 циклов.)

Ответ 7

Вот отличный сайт со многими вариантами: https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/arcsin/onlyelem.html

Лично я пошел в частное приближение Чебышева-Паде со следующим кодом:

double arccos(double x) {
const double pi = 3.141592653;
    return pi / 2 - (.5689111419 - .2644381021*x - .4212611542*(2*x - 1)*(2*x - 1)
         + .1475622352*(2*x - 1)*(2*x - 1)*(2*x - 1))
         / (2.006022274 - 2.343685222*x + .3316406750*(2*x - 1)*(2*x - 1) +
             .02607135626*(2*x - 1)*(2*x - 1)*(2*x - 1));
}

Ответ 8

Быстрая реализация arccosine с точностью до 0,5 градуса может основываться на наблюдении, которая для x в [0,1], acos (x) ≈ √ (2 * (1-x)). Дополнительный масштабный коэффициент повышает точность около нуля. Оптимальный коэффициент можно найти простым бинарным поиском. Отрицательные аргументы обрабатываются в соответствии с acos (-x) = π-acos (x).

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

// Approximate acos(a) with relative error < 5.15e-3
// This uses an idea from Robert Harley posting in comp.arch.arithmetic on 1996/07/12
// https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ
float fast_acos (float a)
{
    const float PI = 3.14159265f;
    const float C  = 0.10501094f;
    float r, s, t, u;
    t = (a < 0) ? (-a) : a;  // handle negative arguments
    u = 1.0f - t;
    s = sqrtf (u + u);
    r = C * u * s + s;  // or fmaf (C * u, s, s) if FMA support in hardware
    if (a < 0) r = PI - r;  // handle negative arguments
    return r;
}

float uint_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof(r));
    return r;
}

int main (void)
{
    double maxrelerr = 0.0;
    uint32_t a = 0;
    do {
        float x = uint_as_float (a);
        float r = fast_acos (x);
        double xx = (double)x;
        double res = (double)r;
        double ref = acos (xx);
        double relerr = (res - ref) / ref;
        if (fabs (relerr) > maxrelerr) {
            maxrelerr = fabs (relerr);
            printf ("xx=% 15.8e  res=% 15.8e  ref=% 15.8e  rel.err=% 15.8e\n",
                    xx, res, ref, relerr);
        }
        a++;
    } while (a);
    printf ("maximum relative error = %15.8e\n", maxrelerr);
    return EXIT_SUCCESS;
}

Выход вышеуказанных тестовых лесов должен выглядеть примерно так:

xx= 0.00000000e+000  res= 1.56272149e+000  ref= 1.57079633e+000  rel.err=-5.14060021e-003
xx= 2.98023259e-008  res= 1.56272137e+000  ref= 1.57079630e+000  rel.err=-5.14065723e-003
xx= 8.94069672e-008  res= 1.56272125e+000  ref= 1.57079624e+000  rel.err=-5.14069537e-003
xx=-2.98023259e-008  res= 1.57887137e+000  ref= 1.57079636e+000  rel.err= 5.14071269e-003
xx=-8.94069672e-008  res= 1.57887149e+000  ref= 1.57079642e+000  rel.err= 5.14075044e-003
maximum relative error = 5.14075044e-003

Ответ 9

    float v,x;
    cin>>x;
    v=x;
    x=(x*3.141592654)/180;
    float a,b,c,d,e,f,g,sum=0;
    a=(x*x)/ar[2];
    b=(x*x*x*x)/ar[4];
    c=(x*x*x*x*x*x)/ar[6];
    d=(x*x*x*x*x*x*x*x)/ar[8];
    sum=1-a+b-c+d;
    cout<<"cos("<<v<<") = "<<sum;