我尝试了millianw 的答案,它给了我 4.5 倍的加速,这太棒了。
但是,millianw 链接到的原始文章计算的是正弦,而不是余弦,并且它的计算方式有所不同。(看起来更简单。)
可以预见的是,15 年后那篇文章的 URL ( http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 ) 今天给出了 404,所以我通过archive.org 获取了它我在这里添加它以供后代使用。
不幸的是,即使这篇文章包含多张图片,archive.org 也只存储了前两张。另外,作者的个人资料页面(http://forum.devmaster.net/users/Nick)没有被存储,所以我想我们永远不会知道尼克是谁。
====================================================
快速准确的正弦/余弦
尼克 06 年 4 月
大家好,
在某些情况下,您需要以非常高的性能运行的正弦和余弦的良好近似值。一个示例是实现圆形表面的动态细分,类似于 Quake 3 中的那些。或者实现波动,以防没有可用的顶点着色器 2.0。
标准的 C sinf() 和 cosf() 函数非常慢,并且提供的精度比我们实际需要的要高得多。我们真正想要的是在精度和性能之间提供最佳折衷的近似值。最著名的近似方法是使用约 0 的泰勒级数(也称为麦克劳林级数),对于正弦,它变为:
x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ...
当我们绘制它时,我们得到:taylor.gif。
绿线是真正的正弦,红线是泰勒级数的前四项。这似乎是一个可以接受的近似值,但让我们仔细看看:taylor_zoom.gif。
在 pi/2 之前它表现得非常好,但之后它很快就会偏离。在 pi 处,它的计算结果为 -0.075 而不是 0。将其用于波浪模拟将导致不可接受的抖动运动。
我们可以添加另一个术语,这实际上可以显着减少误差,但这会使公式变得非常冗长。对于 4 项版本,我们已经需要 7 次乘法和 3 次加法。泰勒系列无法为我们提供我们正在寻找的精度和性能。
然而,我们确实注意到我们需要 sine(pi) = 0。我们可以从 taylor_zoom.gif 中看到另一件事:这看起来非常像抛物线!因此,让我们尝试找到与它尽可能匹配的抛物线公式。抛物线的通用公式是 A + B x + C x\^2。所以这给了我们三个自由度。显而易见的选择是我们想要 sine(0) = 0、sine(pi/2) = 1 和 sine(pi) = 0。这给了我们以下三个等式:
A + B 0 + C 0^2 = 0
A + B pi/2 + C (pi/2)^2 = 1
A + B pi + C pi^2 = 0
其中解 A = 0, B = 4/pi, C = -4/pi\^2。所以我们的抛物线近似变为 4/pi x - 4/pi\^2 x\^2。绘制这个我们得到:抛物线.gif。这看起来比 4 项泰勒级数更糟糕,对吧?错误的!最大绝对误差为 0.056。此外,这种近似将使我们得到平滑的波动,并且只需 3 次乘法和 1 次加法即可计算!
不幸的是,它还不是很实用。这就是我们在 [-pi, pi] 范围内得到的:negative.gif。很明显,我们至少想要一个完整的时期。但同样清楚的是,它只是另一条抛物线,反映在原点周围。它的公式是 4/pi x + 4/pi\^2 x\^2。所以简单的(伪C)解决方案是:
if(x > 0)
{
y = 4/pi x - 4/pi^2 x^2;
}
else
{
y = 4/pi x + 4/pi^2 x^2;
}
但是添加分支不是一个好主意。它使代码显着变慢。但是看看这两个部分到底有多相似。根据 x 的符号,减法变为加法。在第一次尝试消除分支时,我们可以使用 x / abs(x) “提取” x 的符号。除法非常昂贵,但看看得到的公式:4/pi x - x / abs(x) 4/pi\^2 x\^2。通过反转除法,我们可以将其简化为非常漂亮和干净的 4/pi x - 4/pi\^2 x abs(x)。因此,只需一个额外的操作,我们就得到了正弦近似值的两半!这是确认结果的公式图表:abs.gif。
现在让我们看看余弦。基本三角学告诉我们余弦(x) = 正弦(pi/2 + x)。就是这样,将 pi/2 添加到 x 吗?不,我们实际上又得到了抛物线中不需要的部分:shift_sine.gif。我们需要做的是当 x > pi/2 时“环绕”。这可以通过减去 2 pi 来完成。所以代码变成了:
x += pi/2;
if(x > pi) // Original x > pi/2
{
x -= 2 * pi; // Wrap: cos(x) = cos(x - 2 pi)
}
y = sine(x);
又是一个分支。为了消除它,我们可以使用二进制逻辑,如下所示:
x -= (x > pi) & (2 * pi);
请注意,这根本不是有效的 C 代码。但它应该澄清这是如何工作的。当 x > pi 为假时,& 操作将右手部分归零,因此减法不执行任何操作,这是完全等价的。我将把它作为练习留给读者为此创建工作 C 代码(或继续阅读)。显然,余弦比正弦需要更多的操作,但似乎没有任何其他方式,而且它仍然非常快。
现在,0.056 的最大误差很好,但显然 4 项泰勒级数的平均误差仍然较小。回想一下我们的正弦的样子:abs.gif。那么,我们不能做些什么来以最小的成本进一步提高精度吗?当然,当前版本已经适用于许多看起来像正弦的东西与真正的正弦一样好的情况。但对于其他情况,这还不够好。
查看图表,您会注意到我们的近似值总是高估实际正弦值,但 0、pi/2 和 pi 除外。所以我们需要的是在不触及这些重要点的情况下“缩小规模”。解决方案是使用平方抛物线,如下所示:squared.gif。请注意它如何保留这些重要点,但它总是低于真正的正弦值。所以我们可以使用两者的加权平均值来获得更好的近似值:
Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2
Q + P = 1。您可以使用绝对误差或相对误差的精确最小化,但我会为您节省数学。最佳权重是绝对误差的 Q = 0.775,P = 0.225 和相对误差的 Q = 0.782,P = 0.218。我会使用前者。结果图是:average.gif。红线去哪儿了?它几乎完全被绿线覆盖,这立即显示了这个近似值的真实程度。最大误差约为 0.001,提高了 50 倍!公式看起来很长,但括号之间的部分与抛物线的值相同,只需计算一次。事实上,实现这种精度提升只需要 2 次额外的乘法和 2 次加法。
为了使它对负 x 也起作用,我们需要第二个 abs() 操作,这不足为奇。正弦的最终 C 代码变为:
float sine(float x)
{
const float B = 4/pi;
const float C = -4/(pi*pi);
float y = B * x + C * x * abs(x);
#ifdef EXTRA_PRECISION
// const float Q = 0.775;
const float P = 0.225;
y = P * (y * abs(y) - y) + y; // Q * y + P * y * abs(y)
#endif
}
所以我们只需要5次乘法和3次加法;如果我们忽略 abs(),仍然比 4 项 Taylor 更快,而且更精确!余弦版本只需要对 x 进行额外的移位和换行操作。
最后但同样重要的是,如果我不包括 SIMD 优化的汇编版本,我就不会是 Nick。它允许非常有效地执行换行操作,所以我会给你余弦:
// cos(x) = sin(x + pi/2)
addps xmm0, PI_2
movaps xmm1, xmm0
cmpnltps xmm1, PI
andps xmm1, PIx2
subps xmm0, xmm1
// Parabola
movaps xmm1, xmm0
andps xmm1, abs
mulps xmm1, xmm0
mulps xmm0, B
mulps xmm1, C
addps xmm0, xmm1
// Extra precision
movaps xmm1, xmm0
andps xmm1, abs
mulps xmm1, xmm0
subps xmm1, xmm0
mulps xmm1, P
addps xmm0, xmm1
此代码并行计算四个余弦,对于大多数 CPU 架构,每个余弦的峰值性能约为 9 个时钟周期。理想情况下,正弦波只需要 6 个时钟周期。较低精度的版本甚至可以在每个正弦波运行 3 个时钟周期......而且不要忘记 -pi 和 pi 之间的所有输入都是有效的,并且公式在 -pi、-pi/2、0、pi/2 和圆周率。
因此,结论是永远不要再使用泰勒级数来逼近正弦或余弦!为了在本文中添加有用的讨论,我很想知道是否有人知道其他超越函数(如指数、对数和幂函数)的良好近似值。
干杯,
缺口
====================================================
通过访问 Web 存档页面,您可能还会发现本文后面的评论很有趣:
http://web.archive.org/web/20141220225551/http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648