Алгоритм Fast Arc Cos?

Question

У меня есть собственная, очень быстрая функция 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) Спасибо.

Answer 1

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

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

Максимальная погрешность составляет около 0,18 рад.

Answer 2

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

Answer 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

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

Answer 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 1.57079-1.57079*x
end

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

Answer 5

Вы можете аппроксимировать обратный косинус с помощью полинома , как предлагает 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 циклов.)

Answer 6

Другой подход, который вы можете использовать, это использовать комплексные числа.Из формулы Моавра ,

ⅈ ^x = cos (π / 2 * x) + ⅈ * sin (π / 2 * x)

Пусть θ = π / 2 * x.Тогда x = 2θ / π, поэтому

• sin (θ) = ℑ (ⅈ ^ ^{2θ / π} )
• cos (θ) = ℜ (ⅈ^ ^{2θ / π} )

Как вычислить силы powers без греха и кос?Начните с предварительно вычисленной таблицы для степеней 2:

• ⅈ ⁴ = 1
• ⅈ ² = -1
• ⅈ ¹ = ⅈ
• ⅈ ^1/2 = 0,7071067811865476 + 0,7071067811865475 * ⅈ
• ⅈ ^1/4 = 0,9238795325112867 + 0,3826834323650898 * ⅈ
• ⅈ ^1/8 = 0,9807852804032304 + 0,19509032201612825 * ⅈ
• ⅈ ^{1/16 * 10492 0,06 219 096 599}
• ⅈ ^1/32 = 0.9987954562051724 + 0.049067674327418015 * ⅈ
• ⅈ ^1/64 = 0.9996988186962042 + 0.024541228522912288 * 10 *1053*10 ^1/128 = 0.9999247018391445 + 0.012271538285719925 * ⅈ
• ⅈ ^1/256 = 0.9999811752826011 + 0.006135884649154475 * ⅈ

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

Например, tо найти грех и cos 72 ° = 0,8π / 2:

ⅈ ^0,8 ≈ ⅈ ^205/256 = ⅈ ^0b11001101 = ⅈ ^1/2 * ⅈ ^1/4 * ⅈ ^1/32 * ⅈ ^1/64 * ⅈ ^1/256
= 0,3078496400415349 + 0,9514350209690084 * ⅈ

• sin (72 °) ≈ 0,9514350209690084 («точное» значение 0,9510565162951535)
• cos 0 721549 (723064)"точное" значение равно 0,30901699437494745).

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

Например, чтобы найти асин (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 * ⅈ.Грех слишком велик, поэтому уменьшите х на 1/32.
• ⅈ ^13/32 = 0.8032075314806449 + 0.5956993044924334 * ⅈ.

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

Если мы остановимся здесь, мы получим acos (0,6) ≈ 13/32 * π / 2 = 0,6381360077604268 («точное» значение равно 0,6435011087932844.)

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

Answer 7

Быстрая реализация арккозина с точностью около 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's 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

Answer 8

Вот отличный сайт с большим количеством опций: 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));
}