Есть ли точное приближение функции acos()?
Мне нужен acos() Функция с двойной точностью в вычислительном шейдере. Поскольку нет встроенной функции acos() в GLSL с двойной точностью я пытался реализовать свой.
Сначала я реализовал ряд Тейлора, подобный уравнению из ряда Вики - Тейлора, с предварительно вычисленными значениями факультета. Но это кажется неточным в районе 1. Максимальная ошибка была около 0,08 с 40 итерациями.
Я также реализовал этот метод, который очень хорошо работает на процессоре с максимальной ошибкой -2.22045e-16, но у меня есть некоторые проблемы с реализацией этого в шейдере.
В настоящее время я использую acos() функция приближения отсюда, где кто-то разместил свои функции приближения на этом сайте. Я использую самую точную функцию этого сайта, и теперь я получаю максимальную ошибку -7.60454e-08, но эта ошибка слишком высока.
Мой код этой функции:
double myACOS(double x)
{
double part[4];
part[0] = 32768.0/2835.0*sqrt(2.0-sqrt(2.0+sqrt(2.0+sqrt(2.0+2.0*x))));
part[1] = 256.0/135.0*sqrt(2.0-sqrt(2.0+sqrt(2.0+2.0*x)));
part[2] = 8.0/135.0*sqrt(2.0-sqrt(2.0+2.0*x));
part[3] = 1.0/2835.0*sqrt(2.0-2.0*x);
return (part[0]-part[1]+part[2]-part[3]);
}
Кто-нибудь знает другой метод реализации acos() что очень точно и, если возможно, легко реализовать в шейдере?
Некоторая системная информация:
- Nvidia GT 555M
- работает OpenGL 4.3 с optirun
4 ответа
Моя текущая точная реализация шейдера 'acos()' представляет собой смесь из обычной серии Тейлора и ответа Бенса. С 40 итерациями я получаю точность 4.44089e-16 для реализации 'acos()' из math.h. Может быть, это не самый лучший, но у меня работает
И вот оно:
double myASIN2(double x)
{
double sum, tempExp;
tempExp = x;
double factor = 1.0;
double divisor = 1.0;
sum = x;
for(int i = 0; i < 40; i++)
{
tempExp *= x*x;
divisor += 2.0;
factor *= (2.0*double(i) + 1.0)/((double(i)+1.0)*2.0);
sum += factor*tempExp/divisor;
}
return sum;
}
double myASIN(double x)
{
if(abs(x) <= 0.71)
return myASIN2(x);
else if( x > 0)
return (PI/2.0-myASIN2(sqrt(1.0-(x*x))));
else //x < 0 or x is NaN
return (myASIN2(sqrt(1.0-(x*x)))-PI/2.0);
}
double myACOS(double x)
{
return (PI/2.0 - myASIN(x));
}
Любые комментарии, что можно сделать лучше? Например, используя LUT для значений фактора, но в моем шейдере acos () вызывается только один раз, поэтому в этом нет необходимости.
Графический процессор NVIDIA GT 555M - это устройство с вычислительной способностью 2.1, поэтому имеется встроенная аппаратная поддержка для базовых операций с двойной точностью, в том числе с плавным сложением (FMA). Как и во всех графических процессорах NVIDIA, эмуляция квадратного корня. Я знаком с CUDA, но не GLSL. Согласно версии 4.3 спецификации GLSL, она представляет FMA двойной точности как функцию fma() и обеспечивает квадратный корень двойной точности, sqrt(), Не ясно, является ли sqrt() реализация правильно округлена в соответствии с правилами IEEE-754. Я предполагаю, что это по аналогии с CUDA.
Вместо использования ряда Тейлора можно было бы использовать полиномиальное минимаксное приближение, тем самым уменьшая количество требуемых членов. Минимаксные приближения обычно генерируются с использованием варианта алгоритма Ремеза. Для оптимизации как скорости, так и точности, использование FMA имеет важное значение. Оценка полинома по схеме Хорнера способствует высокой степени накопления. В приведенном ниже коде используется схема Хорнера второго порядка. Как и в ответе DanceIgel, acosудобно рассчитывается с использованиемasinаппроксимация как основной строительный блок в сочетании со стандартными математическими тождествами.
Для 400M векторов испытаний максимальная относительная ошибка, наблюдаемая с помощью приведенного ниже кода, составляла 2.67e-16, а максимальная наблюдаемая ошибка ulp - 1.442 ulp.
/* compute arcsin (a) for a in [-9/16, 9/16] */
double asin_core (double a)
{
double q, r, s, t;
s = a * a;
q = s * s;
r = 5.5579749017470502e-2;
t = -6.2027913464120114e-2;
r = fma (r, q, 5.4224464349245036e-2);
t = fma (t, q, -1.1326992890324464e-2);
r = fma (r, q, 1.5268872539397656e-2);
t = fma (t, q, 1.0493798473372081e-2);
r = fma (r, q, 1.4106045900607047e-2);
t = fma (t, q, 1.7339776384962050e-2);
r = fma (r, q, 2.2372961589651054e-2);
t = fma (t, q, 3.0381912707941005e-2);
r = fma (r, q, 4.4642857881094775e-2);
t = fma (t, q, 7.4999999991367292e-2);
r = fma (r, s, t);
r = fma (r, s, 1.6666666666670193e-1);
t = a * s;
r = fma (r, t, a);
return r;
}
/* Compute arccosine (a), maximum error observed: 1.4316 ulp
Double-precision factorization of π courtesy of Tor Myklebust
*/
double my_acos (double a)
{
double r;
r = (a > 0.0) ? -a : a; // avoid modifying the "sign" of NaNs
if (r > -0.5625) {
/* arccos(x) = pi/2 - arcsin(x) */
r = fma (9.3282184640716537e-1, 1.6839188885261840e+0, asin_core (r));
} else {
/* arccos(x) = 2 * arcsin (sqrt ((1-x) / 2)) */
r = 2.0 * asin_core (sqrt (fma (0.5, r, 0.5)));
}
if (!(a > 0.0) && (a >= -1.0)) { // avoid modifying the "sign" of NaNs
/* arccos (-x) = pi - arccos(x) */
r = fma (1.8656436928143307e+0, 1.6839188885261840e+0, -r);
}
return r;
}
«Проблема» в том, что у него есть особенности (точки ветвления) при ±1, и по этой причине Тейлор или Паде не могут хорошо аппроксимироваться. То же самое относится и к приближениям MiniMax, которые пытаются минимизировать максимальную норму, когда область содержит особенность.
Сингулярности в основном представляют собой квадратные корни, поэтому их можно исключить, а оставшаяся функция будет гладкой и хорошо аппроксимирует. Оставшаяся функция находится в следующем коде C99, это рациональная аппроксимация MiniMax, вычисленная с помощью алгоритма Ремеза . (На самом деле код GNU-C99, поэтомуM_PIдоступен).
Примечания к коду:
Аппроксимация минимизирует относительную ошибку, поскольку речь идет о числах с плавающей запятой.
Он использует (квадратный корень),
fma(плавающее умножение-сложение),ldexp(масштабировать в степени 2)Оно использует
if-elseи условное присвоение. Я не знаком с шейдерами и не влияет ли это на производительность. Если это так, возможно, код можно реорганизовать. Каждый путь развиваетasinQдля некоторого значения 0 ≤ x ≤ 1/2.
#include <math.h>
#include <stdio.h>
static const double P[] =
{
-5254.7508920534192751630665,
+6236.4128689522053163220058,
-2597.2384260994230261835544,
+445.05537923135581032176275,
-27.342673770163272135013850,
+0.30325426599090297962612977
};
static const double Q[] =
{
-5254.7508920534192747096309,
+6674.3087766233233966852037,
-3054.9042449253514128522165,
+603.81083008747677276795370,
-47.647718147243950851054008
// 1.0
};
// Approximate arcsin(sqrt(x/2)) / sqrt(x/2) as MiniMax [5/5] over [0, 1/2].
static inline double asinQ (double x)
{
double p = fma (P[5], x, P[4]);
p = fma (p, x, P[3]);
p = fma (p, x, P[2]);
p = fma (p, x, P[1]);
p = fma (p, x, P[0]);
double q = Q[4] + x; // Q[5] = 1.0
q = fma (q, x, Q[3]);
q = fma (q, x, Q[2]);
q = fma (q, x, Q[1]);
q = fma (q, x, Q[0]);
return p / q;
}
double my_acos (double x)
{
double x_abs = x > 0.0 ? x : -x;
if (x_abs <= 0.5)
{
return M_PI/2 - x * asinQ (ldexp (x*x, 1));
}
else
{
double x1 = 1.0 - x_abs;
double ret = sqrt (ldexp (x1, 1)) * asinQ (x1);
return x > 0.0 ? ret : M_PI - ret;
}
}
int main (int argc, char *argv[])
{
double max_err = -1.0;
double min_bits = 0;
double max_x = 0;
int ex = 8;
if (argc >= 2)
sscanf (argv[1], "%i", &ex);
for (double x = -1.0; x <= 1.0; x += ldexp (1, -ex))
{
double err = (my_acos (x) - acos(x)) / acos (x);
double bits = - log2 (fabs (err));
//printf ("%+.6f: % 6e (%.2f bits)\n", x, err, bits);
if (fabs (err) > max_err)
max_err = fabs (err), min_bits = bits, max_x = x;
}
printf ("worst (%ld values): x=%+.6f, err=%6e, bits=%.2f\n",
1 + (2L << ex), max_x, max_err, min_bits);
return 0;
}
Запуск этого для отпечатков значений 16 млн.
worst (16777217 values): x=+0.839781, err=5.803408e-16, bits=50.61
Таким образом, максимальная относительная ошибка составляет около 6·10−16 , что соответствует точности около 50,5 бит. (Абсолютная ошибка для этого значения равна 3,4·10 −16 .) Полагаяacosиsqrtввести ошибку 0,5 ULP, точностьmy_acosсоставляет около 1,5 ULP.
Возможно, это решение поможет. Это лучше, чем 1% правильного угла для x = 1 до 0,2.
acos(x) ~= sqrt(2) * ( sqrt(1-x) + (1/11) * (sqrt(1-x))^3 )
Это началось с серии Тейлора, предоставленной Wolfram. Это требовало слишком много терминов даже для приблизительного значения ниже 0,8. В этом методе использовалась общая форма первых двух членов, но были изменены коэффициенты для улучшения соответствия. Интересно, что сработал целочисленный коэфф целого числа 11.