三角函数逼近快速算法(正余弦).doc

三角函数逼近快速算法(正余弦) 原文出自： http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/ 里面提到了查表，采用查表并配合插值；以及泰勒级数 看过第一篇的文章后，大呼过瘾！原文作者的思路非常简捷，有趣，偶英语比较差，欢迎指正，废话不多说看文章 原文出处： /forums/showthread.php?t=5784 http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/ 在某些情况下我们需要一些更高效的且近似于标准值的 sin 和cos函数。 有时候我们并需要过高的精度，这时 C语言中自带的三角函数（sinf() 和 cosf() f）计算的精度超出了我们所需要的精度要求，所以其效率很低。我们真正需要的是间于精度和效率的一个折中的方案。众所周知的取近似值的方法是：泰勒级数（和著名的马克劳林级数） 代码是： x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ... 我们绘制了下图: 绿线是标准的sin函数，红线是4项泰勒级数展开式。这个近似值的效果看起来还不错，但是如果你仔细观察后会发现 它在pi/2之前的效果还是很好的，但是超过了pi/2后就开始快速偏离标准sin。它在pi处的值比原来的0多了0.075.用这个方法来模拟波动忽动忽停，看起来很呆板，这个效果当然不是我们想要的。 我们可以继续增加泰勒级数项的个数来减小误差，但是这将导致我们的公式非常的冗长。用4项的泰勒级数展开式需要我们进行7次乘法和3次加法来完成。泰勒级数不能同时满足我对精度和效率的要求。 刚刚近似如果能满足sin（pi）=0就好了。从上图我还可以发现另一件事：这个曲线看起来很象抛物线！所以我们来寻找一个尽可能和sin接近的抛物线（公式）。抛物线的范式方程是：A + B x + C x^2.这个公式可以让们控制三个自由度。显然我们需要其满足sine(0) = 0, sine(pi/2) = 1 and sine(pi) = 0. 这样我们就得到了3个等式。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.我们的抛物线诞生啦！ 貌似这个的误差看起来比泰勒级数还要遭。其实不是的！这种方法的最大误差是0.056.（译者：而且这种近似值没有误差积累）而且这个近似值的绘制出的波动是光滑的，而且只需要3次乘法和一次加法。不过它还不够完美。下图是[-pi, pi] 之间的图像： 显然我们至少需要它在1个完整的周期内都符合我们要求。但是我们可以看出，我们缺少的另一半是原抛物线的一个映射。它的公式是：4/pi x + 4/pi^2 x^2。所以我们可以直接这样写： Code: 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)。除法的代价是非常大的，我们来观察一下现在的公式： 4/pi x - x / abs(x) 4/pi^2 x^2。将除法化简后我们得到：4/pi x - 4/pi^2 x abs(x).所以只需要多一步运算我们就得到了我们需要的sin逼近值！下图是结果 接下来我们要考虑cos。有基础的三角公理可以知道：cos(x) = sin(pi/2 + x).把x多加一个pi/2就可以搞定了？事实上它的某一部分不是我们期望得到的。 我们需要做的就是当x pi/2时“跳回”。这个可以由减去2 pi来实现。 Code: 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); Code: x -= (x pi) (2 * pi); 注意这并不是c的源代码。但是这个应该可以说明它是怎么样运行的。当x pi是false 时，逻辑“与”（）运算后得到的