acos() 函数是否有准确的近似值？

Question

我需要一个在计算着色器中具有双精度的acos() 函数。由于GLSL中没有双精度的acos()的内置函数，我尝试自己实现。

起初，我实现了一个泰勒级数，例如 Wiki - Taylor series 中的方程，并带有预先计算的教师值。但这似乎在 1 左右不准确。最大误差约为 0.08，迭代 40 次。

我还实现了this method，它在 CPU 上运行良好，最大误差为 -2.22045e-16，但在着色器中实现它时遇到了一些麻烦。

目前，我正在使用来自here 的acos() 近似函数，其中有人在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() 的另一种实现方法吗？它非常准确并且——如果可能的话——很容易在着色器中实现？

一些系统信息：

• 英伟达 GT 555M
• 使用 optirun 运行 OpenGL 4.3

Answer 1

NVIDIA GT 555M GPU 是一款具有 2.1 计算能力的设备，因此对基本的双精度运算有原生硬件支持，包括fused multipy-add (FMA)。与所有 NVIDIA GPU 一样，平方根运算是模拟的。我熟悉 CUDA，但不熟悉 GLSL。根据GLSL specification 的4.3 版，它将双精度FMA 公开为函数fma() 并提供双精度平方根sqrt()。目前尚不清楚sqrt() 实现是否根据IEEE-754 规则正确舍入。我会假设它是，与 CUDA 类比。

我们希望使用多项式minimax approximation，而不是使用泰勒级数，从而减少所需的项数。 Minimax 近似值通常使用Remez algorithm 的变体生成。为了优化速度和准确性，使用 FMA 是必不可少的。使用Horner scheme 对多项式进行评估有助于提高准确性。在下面的代码中，使用了二阶霍纳方案。正如在 DanceIgel 的 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;
}

Answer 2

我当前对“acos()”的准确着色器实现是通常的泰勒级数和Bence 的答案的混合。通过 40 次迭代，我从 math.h 中获得了 'acos()' 实现的 4.44089e-16 精度。也许它不是最好的，但它对我有用：

这里是：

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));
}

任何cmets，有什么可以做得更好？例如，对因子的值使用 LUT，但在我的着色器中，'acos()' 只调用一次，因此不需要它。