Есть ли точная аппроксимация функции acos()?

Question

Мне нужна функция acos() с двойной точностью в вычислительном шейдере. Поскольку встроенная функция acos() в GLSL с двойной точностью, я попытался реализовать свои собственные.

Сначала я реализовал серию Тейлора, как уравнение от Wiki - Taylor series с заранее рассчитанными значениями факультета. Но это кажется неточным около 1. Максимальная ошибка была примерно 0,08 с 40 итерациями.

Я также реализовал this method, который очень хорошо работает на процессоре с максимальной ошибкой -2.22045e-16, но у меня есть некоторые проблемы для реализации этого в шейдере.

В настоящее время я использую функцию приближения acos() от here, где кто-то разместил свои функции аппроксимации на сайте this. Я использую самую точную функцию этого сайта, и теперь я получаю максимальную ошибку -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(), который является очень точным и -если possible- легко реализовать в шейдере?

Определенная информация о системе:

• Nvidia GT 555M
• работает OpenGL 4.3 с optirun

Answer 1

Моя текущая точная реализация шейдера из 'экоса()' представляет собой смесь из обычного Тейлором серия и ответ от Bence. С 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() просто вызывается один раз, поэтому нет необходимости в нем.

Answer 2

Графический процессор NVIDIA GT 555M - это устройство с вычислительной способностью 2.1, поэтому имеется базовая аппаратная поддержка основных операций с двойной точностью, включая fused multipy-add (FMA). Как и во всех графических процессорах NVIDIA, эскалация квадратного корня. Я знаком с CUDA, но не с GLSL. Согласно версии 4.3 GLSL specification, он предоставляет функцию FMA с двойной точностью как функцию fma() и обеспечивает квадратный корень с двойной точностью, sqrt(). Неясно, правильно ли выполнена реализация sqrt() в соответствии с правилами IEEE-754. Я предполагаю, что это, по аналогии с CUDA.

Вместо использования серии Тейлора можно было бы использовать многочлен minimax approximation, тем самым уменьшив количество требуемых терминов. Минимаксные аппроксимации обычно генерируются с использованием варианта Remez algorithm. Для оптимизации скорости и точности использование FMA имеет важное значение. Оценка многочлена с Horner scheme является достаточной для высокого начисления. В приведенном ниже коде используется схема Орнера второго порядка.Как и в файле DanceIgel's answer, 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; 
}