Есть ли точное приближение функции acos()?

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

Вместо использования ряда Тейлора можно было бы использовать полиномиальное минимаксное приближение, тем самым уменьшая необходимое количество терминов. Минимаксные приближения обычно генерируются с использованием варианта алгоритма Ремеза. Для оптимизации скорости и точности необходимо использовать FMA. Оценка полинома с помощью схемы Хорнера способствует высокой точности. В приведенном ниже коде используется схема Хорнера второго порядка. Как и в ответе от DanceIgel, acos удобно вычислять с использованием приближения asin в качестве базового строительного блока в сочетании со стандартными математическими тождествами.

При использовании 400 млн тестовых векторов максимальная относительная ошибка, наблюдаемая в приведенном ниже коде, составила 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;
}

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