我需要在计算着色器中具有双精度的
acos()acos()首先，我使用预先计算的教师值实现了一个泰勒级数，如 Wiki - 泰勒级数 中的方程。但这似乎在 1 左右不准确。迭代 40 次时，最大误差约为 0.08。
我还实现了这个方法，它在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()一些系统信息：
NVIDIA GT 555M GPU 是计算能力为 2.1 的设备，因此对基本双精度运算有原生硬件支持，包括 fused multipy-add (FMA)。与所有 NVIDIA GPU 一样，平方根运算是模拟的。我熟悉CUDA，但不熟悉GLSL。根据 GLSL 规范 4.3 版，它将双精度 FMA 作为函数
fma()sqrt()sqrt()不使用泰勒级数，而是使用多项式极大极大近似，从而减少所需的项数。极小极大近似通常是使用 Remez 算法 的变体生成的。为了优化速度和准确性，FMA 的使用至关重要。使用Horner 方案 评估多项式有助于获得较高的准确度。在下面的代码中，使用了二阶霍纳方案。正如 DanceIgel 的 answer 中一样，
acosasin使用 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;
}
我当前的“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));
}
有什么意见，可以做得更好吗？例如，使用 LUT 作为因子的值，但在我的着色器中“acos()”仅被调用一次，因此不需要它。
也许这个解决方案有帮助。对于 x = 1 到 0.2，它比正确角度的 1% 更好。
acos(x) ~= sqrt(2) * ( sqrt(1-x) + (1/11) * (sqrt(1-x))^3 )
这始于 Wolfram 提供的泰勒级数。即使粗略值低于 0.8，也需要太多项。该方法使用前两项的一般形式，但更改了系数以改进匹配。有趣的是，整数 11 的整数系数有效。
acos奇点基本上都是平方根，因此可以将它们分解出来，剩下的函数是平滑的并且可以很好地近似。剩下的函数是后面的 C99 代码中的
asinQM_PI代码注释：
近似值最小化了相对误差，因为它是关于浮点的。
它使用
sqrtfmaldexp它使用
if-elseasinQ#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 个 Mio 值
worst (16777217 values): x=+0.839781, err=5.803408e-16, bits=50.61
因此最大相对误差约为 6·10−16，精度约为 50.5 位。 （该值的绝对误差为 3.4·10−16。）假设
acossqrtmy_acos