快速弧算法算法？
jma*_*erx 22 c c++ algorithm math performance
我有自己的,非常快的cos函数:
float sine(float x)
{
const float B = 4/pi;
const float C = -4/(pi*pi);
float y = B * x + C * x * abs(x);
// const float Q = 0.775;
const float P = 0.225;
y = P * (y * abs(y) - y) + y; // Q * y + P * y * abs(y)
return y;
}
float cosine(float x)
{
return sine(x + (pi / 2));
}
但是现在当我描述时,我发现acos()正在杀死处理器.我不需要高度精确.什么是计算acos的快速方法(x)谢谢.
dan*_*n04 39
一个简单的三次近似,x∈{-1,-½,0,½,1}的拉格朗日多项式是:
double acos(x) {
return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966;
}
它的最大误差约为0.18rad.
- 最大误差为 10.31 度。相当大，但在某些解决方案中可能就足够了。适用于计算速度比精度更重要的情况。四次近似可能会产生更高的精度并且仍然比本地 acos 更快吗？ (2认同)
spe*_*der 22
有多余的记忆？查找表(如果需要,带插值)将是最快的.
- @Jex:bounds - 检查你的参数(它必须在-1和1之间).然后乘以2的漂亮幂,比如64,产生范围(-64,64).添加64使其为非负(0,128).使用整数部分索引查找表,如果需要,使用小数部分在两个最接近的条目之间进行插值.如果您不想插值,请尝试添加64.5并发言,这与舍入到最近相同. (7认同)
- 查找表需要一个索引,这将需要一个float到int的转换,这可能会破坏性能. (4认同)
- @phkahler：float 到 int 的转换在 x86 上非常便宜，几乎与 FP 添加一样便宜，正如您在 Agner Fog 的延迟/吞吐量/uop 表中看到的那样（http://agner.org/optimize）。对索引进行范围检查以确保它不会在表之外建立索引可能同样昂贵。`int idx = x * 4096.0` 在 Intel Haswell 上将有大约 9 个周期的延迟。到目前为止，其中最昂贵的部分是来自适当大小的表的缓存未命中。如果没有一堆不依赖于 acos 结果的并行计算，大表可能会更慢（尤其是在缓存竞争的情况下）。 (2认同)
Fno*_*ord 16
nVidia有一些很好的资源,展示了如何近似非常昂贵的数学函数,例如:acos asin atan2 等等......
当执行速度比精度更重要(在合理范围内)时,这些算法会产生良好的结果.这是他们的acos功能:
// Absolute error <= 6.7e-5
float acos(float x) {
float negate = float(x < 0);
x = abs(x);
float ret = -0.0187293;
ret = ret * x;
ret = ret + 0.0742610;
ret = ret * x;
ret = ret - 0.2121144;
ret = ret * x;
ret = ret + 1.5707288;
ret = ret * sqrt(1.0-x);
ret = ret - 2 * negate * ret;
return negate * 3.14159265358979 + ret;
}
以下是计算acos(0.5)时的结果:
nVidia: result: 1.0471513828611643
math.h: result: 1.0471975511965976
那非常接近!根据您所需的精确度,这可能是一个很好的选择.
- 与 nvidia 网站上的“参考实现”中的评论相反，绝对错误不 <= 6.7e-5，但我能够观察到 6.759167e-05 的错误。另外，你有多确定这个函数实际上比普通的 acos 更快？在第四代核心 i5 (Haswell) 上，nvidia 函数始终比 math.h 中的 acos 慢 25%。 (2认同)
- 是否有理由使用`3.14159265358979`代替`math.pi`？ (2认同)
- 提出问题的原因是导致问题的差别很小：将其移植到任何其他语言时-是否应使用pi的常量？还是这是有意修改的pi，调整为与近似值更好地工作？*（可以进行测试，当然不会有太多麻烦）* (2认同)
我有我自己的.它非常准确,速度很快.它取决于我围绕四次收敛建立的定理.它真的很有趣,你可以看到方程式以及它使我的自然对数近似收敛的速度:https://www.desmos.com/calculator/yb04qt8jx4
这是我的arccos代码:
function acos(x)
local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
local c=0.88-0.77*x c=(c+(2-a)/c)/2
return (8*(c+(2-a)/c)-(b+(2-2*x)/b))/6
end
很多只是平方根近似.它的效果也非常好,除非你太接近于取0的平方根.它的平均误差(不包括x = 0.99到1)为0.0003.但问题是,在0.99时,它开始变得很糟糕,而在x = 1时,精度的差异变为0.05.当然,这可以通过在平方根(lol nope)上进行更多迭代来解决,或者只是一点点,如果x> 0.99然后使用一组不同的平方根线性化,但这会使代码变得冗长而丑陋.
如果你不太关心准确性,那么你可以在每个平方根上进行一次迭代,这仍然可以保持在0.0162范围内的某个位置,或者精度如下:
function acos(x)
local a=1.43+0.59*x a=(a+(2+2*x)/a)/2
local b=1.65-1.41*x b=(b+(2-2*x)/b)/2
local c=0.88-0.77*x c=(c+(2-a)/c)/2
return 8/3*c-b/3
end
如果您对它没问题,可以使用预先存在的平方根代码.它将摆脱在x = 1时有点疯狂的等式:
function acos(x)
local a = math.sqrt(2+2*x)
local b = math.sqrt(2-2*x)
local c = math.sqrt(2-a)
return 8/3*d-b/3
end
坦率地说,如果你真的有时间紧迫,请记住你可以将arccos线性化为3.14159-1.57079x并且只做:
function acos(x)
return 1.57079-1.57079*x
end
无论如何,如果你想查看我的arccos近似方程式列表,你可以访问https://www.desmos.com/calculator/tcaty2sv8l我知道我的近似值对于某些事情来说不是最好的,但如果你'如果我的近似值很有用,请使用它们,但请尽量给予我信任.
您可以使用dan04建议的多项式来近似反余弦,但是多项式是-1和1附近非常差的近似,其中反余弦的导数变为无穷大.当你增加多项式的次数时,你会迅速达到递减收益,并且仍然很难在端点周围得到一个很好的近似值.一个有理函数(两个多项式的商)可以在这个情况下,一个更好的近似.
acos(x) ? ?/2 + (ax + bx³) / (1 + cx² + dx?)
哪里
a = -0.939115566365855
b = 0.9217841528914573
c = -1.2845906244690837
d = 0.295624144969963174
在区间(-1,1)上的最大绝对误差为0.017弧度(0.96度).下面是一个图(黑色的反余弦,红色的三次多项式近似,上面的蓝色函数)用于比较:
选择上述系数以最小化整个域上的最大绝对误差.如果您愿意在端点处允许更大的错误,则可以使间隔(-0.98,0.98)上的错误小得多.度数为5的分子和度数为2的分母大约与上述函数一样快,但稍微不准确.以性能为代价,您可以通过使用更高次多项式来提高精度.
关于性能的注意事项:计算两个多项式仍然非常便宜,您可以使用融合乘法加法指令.划分并不是那么糟糕,因为你可以使用硬件倒数近似和乘法.与acos近似中的误差相比,倒数近似中的误差可以忽略不计.在2.6 GHz Skylake i7上,这种近似可以使用AVX每6个周期执行大约8个反余弦.(即吞吐量,延迟超过6个周期.)
您可以采用的另一种方法是使用复数。根据de Moivre的公式，
？x = cos（π/ 2 * x）+？* sin（π/ 2 * x）
令θ=π/ 2 * x。那么x =2θ/π
- sin（θ）=？（？^ 2θ/π）
- cos（θ）=？（？^ 2θ/π）
您如何计算的幂？没有罪恶和COS？从一个预先计算的表开始，得到2的幂：
- ？4 = 1
- ？2 = -1
- ？1 =？
- ？1/2 = 0.7071067811865476 + 0.7071067811865475 *？
- ？1/4 = 0.9238795325112867 + 0.3826834323650898 *？
- ？1/8 = 0.9807852804032304 + 0.19509032201612825 *？
- ？1/16 = 0.9951847266721969 + 0.0980171403295606 *？
- ？1/32 = 0.9987954562051724 + 0.049067674327418015 *？
- ？1/64 = 0.9996988186962042 + 0.024541228522912288 *？
- ？1/128 = 0.9999247018391445 + 0.012271538285719925 *？
- ？1/256 = 0.9999811752826011 + 0.006135884649154475 *？
计算？的任意值 x，将指数近似为二进制分数，然后将表中的相应值相乘。
例如，要找到72°= 0.8？/ 2的正弦和余弦：
？0.8
≈？205/256
=？0b11001101
=？1/2 *？1/4 *？1/32 *？1/64 *？1/256
= 0.3078496400415349 + 0.9514350209690084 *？
- sin（72°）≈0.9514350209690084（“精确”值是0.9510565162951535）
- cos（72°）≈0.3078496400415349（“精确”值是0.30901699437494745）。
要查找asin和acos，可以将此表与Bisection方法一起使用：
例如，要找到asin（0.6）（3-4-5三角形中的最小角度）：
- ？0 = 1 + 0 *？正弦太小，因此将x增加1/2。
- ？1/2 = 0.7071067811865476 + 0.7071067811865475 *？。正弦太大，因此将x减小1/4。
- ？1/4 = 0.9238795325112867 + 0.3826834323650898 *？正弦太小，因此将x增加1/8。
- ？3/8 = 0.8314696123025452 + 0.5555702330196022 *？正弦仍然太小，因此将x增加1/16。
- ？7/16 = 0.773010453362737 + 0.6343932841636455 *？正弦太大，因此将x减小1/32。
- ？13/32 = 0.8032075314806449 + 0.5956993044924334 * ?.
每次增加x时，乘以相应的？幂。每次减少x，除以相应的？幂。
如果在此处停止，我们将获得acos（0.6）≈13/32 *π/ 2 = 0.6381360077604268（“精确”值为0.6435011087932844。）
当然，准确性取决于迭代次数。对于快速的近似值，请使用10次迭代。对于“高精确度”，请使用50-60次迭代。
快速反余弦实现，精确到约 0.5 度，可以基于 观察[0,1] 中的 x，acos(x) \xe2\x89\x88 \xe2\x88\x9a(2*(1-x ））。额外的比例因子可提高接近于零的精度。可以通过简单的二分搜索找到最佳因子。负参数按照 acos (-x) = \xcf\x80 - acos (x) 处理。
\n\n#include <stdio.h>\n#include <stdlib.h>\n#include <stdint.h>\n#include <string.h>\n#include <math.h>\n\n// Approximate acos(a) with relative error < 5.15e-3\n// This uses an idea from Robert Harley\'s posting in comp.arch.arithmetic on 1996/07/12\n// https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ\nfloat fast_acos (float a)\n{\n const float PI = 3.14159265f;\n const float C = 0.10501094f;\n float r, s, t, u;\n t = (a < 0) ? (-a) : a; // handle negative arguments\n u = 1.0f - t;\n s = sqrtf (u + u);\n r = C * u * s + s; // or fmaf (C * u, s, s) if FMA support in hardware\n if (a < 0) r = PI - r; // handle negative arguments\n return r;\n}\n\nfloat uint_as_float (uint32_t a)\n{\n float r;\n memcpy (&r, &a, sizeof(r));\n return r;\n}\n\nint main (void)\n{\n double maxrelerr = 0.0;\n uint32_t a = 0;\n do {\n float x = uint_as_float (a);\n float r = fast_acos (x);\n double xx = (double)x;\n double res = (double)r;\n double ref = acos (xx);\n double relerr = (res - ref) / ref;\n if (fabs (relerr) > maxrelerr) {\n maxrelerr = fabs (relerr);\n printf ("xx=% 15.8e res=% 15.8e ref=% 15.8e rel.err=% 15.8e\\n",\n xx, res, ref, relerr);\n }\n a++;\n } while (a);\n printf ("maximum relative error = %15.8e\\n", maxrelerr);\n return EXIT_SUCCESS;\n}\n上述测试脚手架的输出应类似于以下内容：
\n\nxx= 0.00000000e+000 res= 1.56272149e+000 ref= 1.57079633e+000 rel.err=-5.14060021e-003\nxx= 2.98023259e-008 res= 1.56272137e+000 ref= 1.57079630e+000 rel.err=-5.14065723e-003\nxx= 8.94069672e-008 res= 1.56272125e+000 ref= 1.57079624e+000 rel.err=-5.14069537e-003\nxx=-2.98023259e-008 res= 1.57887137e+000 ref= 1.57079636e+000 rel.err= 5.14071269e-003\nxx=-8.94069672e-008 res= 1.57887149e+000 ref= 1.57079642e+000 rel.err= 5.14075044e-003\nmaximum relative error = 5.14075044e-003\n