У меня есть моя, очень быстрая функция 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) Спасибо.
Простая кубическая аппроксимация, полином Лагранжа для x ∈ {-1, -½, 0, ½, 1} является:
double acos(x) {
return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966;
}
Он имеет максимальную погрешность около 0,18 рад.
Есть запасная память? Таблица поиска (с интерполяцией, если требуется) будет самой быстрой.
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
Это довольно близко! В зависимости от вашей требуемой степени точности это может быть хорошим вариантом для вас.
У меня есть своя. Это довольно точно и быстро. Это работает с одной теоремой, построенной вокруг квадрической сходимости. Это действительно интересно, и вы можете видеть уравнение и как быстро это может привести к схождению моего естественного логарифмического приближения: 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. Я знаю, что мои приближения не самые лучшие для некоторые вещи, но если вы делаете что-то, где мои аппроксимации будут полезны, используйте их, но постарайтесь дать мне кредит.
Другой подход, который вы можете предпринять, - использовать сложные числа. Из de Moivre formula,
ⅈ x= cos (& pi;/2 * x) + ⅈ * sin (& pi;/2 * x)
Пусть & theta; = & pi;/2 * x. Тогда x = 2 & theta;/& pi, поэтому
Как вы можете вычислить степени ⅈ без sin и cos? Начните с предварительно вычисленной таблицы для степеней 2:
Чтобы вычислить произвольные значения ⅈ 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 * ⅈ
Чтобы найти asin и acos, вы можете использовать эту таблицу с помощью метода Bisection:
Например, чтобы найти asin (0.6) (наименьший угол в 3-4-5 треугольниках):
Каждый раз, когда вы увеличиваете x, умножьте на соответствующую мощность ⅈ. Каждый раз, когда вы уменьшаете x, делите на соответствующую мощность ⅈ.
Если мы остановимся здесь, получим acos (0.6) & approx; 13/32 * & pi;/2 = 0,6381360077604268 ( "Точное" значение 0,6435011087932844.)
Точность, конечно, зависит от количества итераций. Для быстрого и грязного приближения используйте 10 итераций. Для "интенсивной точности" используйте 50-60 итераций.
Вы можете аппроксимировать обратный косинус с помощью полинома как предложено 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 циклов.)
Вот отличный сайт со многими вариантами: 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));
}
Быстрая реализация 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
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;