opengl - 是否有 acos() 函数的精确近似值？

Question

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

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

我也实现了 this method它在 CPU 上运行良好，最大误差为 -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()的另一种实现方法吗？这是非常准确的 - 如果可能 - 易于在着色器中实现？

一些系统信息:

英伟达 GT 555M

使用 optirun 运行 OpenGL 4.3

Answer 1

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

与其使用泰勒级数，不如使用多项式 minimax approximation ，从而减少所需的术语数量。极小极大近似通常使用 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;
}