《三角函数逼近快速算法(正余弦).doc》由会员分享，可在线阅读，更多相关《三角函数逼近快速算法(正余弦).doc（21页珍藏版）》请在三一办公上搜索。
1、三角函数逼近快速算法(正余弦)原文出自： http:/lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/ 里面提到了查表，采用查表并配合插值；以及泰勒级数 看过第一篇的文章后，大呼过瘾！原文作者的思路非常简捷，有趣，偶英语比较差，欢迎指正，废话不多说看文章 原文出处： http:/lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/ 在某些情况下我们需要一些更高效的且近似于标准值的 sin 和cos函数。 有时候我
2、们并需要过高的精度，这时 C语言中自带的三角函数（sinf() 和 cosf() f）计算的精度超出了我们所需要的精度要求，所以其效率很低。我们真正需要的是间于精度和效率的一个折中的方案。众所周知的取近似值的方法是：泰勒级数（和著名的马克劳林级数） 代码是： x - 1/6 x3 + 1/120 x5 - 1/5040 x7 + . 我们绘制了下图: 绿线是标准的sin函数，红线是4项泰勒级数展开式。这个近似值的效果看起来还不错，但是如果你仔细观察后会发现 它在pi/2之前的效果还是很好的，但是超过了pi/2后就开始快速偏离标准sin。它在pi处的值比原来的0多了0.075.用这个方法来模拟波
3、动忽动忽停，看起来很呆板，这个效果当然不是我们想要的。 我们可以继续增加泰勒级数项的个数来减小误差，但是这将导致我们的公式非常的冗长。用4项的泰勒级数展开式需要我们进行7次乘法和3次加法来完成。泰勒级数不能同时满足我对精度和效率的要求。 刚刚近似如果能满足sin（pi）=0就好了。从上图我还可以发现另一件事：这个曲线看起来很象抛物线！所以我们来寻找一个尽可能和sin接近的抛物线（公式）。抛物线的范式方程是：A + B x + C x2.这个公式可以让们控制三个自由度。显然我们需要其满足sine(0) = 0, sine(pi/2) = 1 and sine(pi) = 0. 这样我们就得到了3
4、个等式。A + B 0 + C 02 = 0 A + B pi/2 + C (pi/2)2 = 1 A + B pi + C pi2 = 0 解得：A = 0, B = 4/pi, C = -4/pi2.我们的抛物线诞生啦！ 貌似这个的误差看起来比泰勒级数还要遭。其实不是的！这种方法的最大误差是0.056.（译者：而且这种近似值没有误差积累）而且这个近似值的绘制出的波动是光滑的，而且只需要3次乘法和一次加法。不过它还不够完美。下图是-pi, pi 之间的图像： 显然我们至少需要它在1个完整的周期内都符合我们要求。但是我们可以看出，我们缺少的另一半是原抛物线的一个映射。它的公式是：4/pi x
5、+ 4/pi2 x2。所以我们可以直接这样写： Code: if(x 0) y = 4/pi x - 4/pi2 x2; else y = 4/pi x + 4/pi2 x2; 添加一个条件分支不是一个好的方法。它会让程序渐渐的变慢。但是观察一下我们模拟的和标准的图像是多么的接近啊！观察上面两式子，只是中间的一项正负号不同，我的第一个单纯的想法是可以提取x的正负号来消除分支，即使用：x / abs(x)。除法的代价是非常大的，我们来观察一下现在的公式： 4/pi x - x / abs(x) 4/pi2 x2。将除法化简后我们得到：4/pi x - 4/pi2 x abs(x).所以只需要多一
6、步运算我们就得到了我们需要的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 * p
7、i); Code: x -= (x pi) & (2 * pi); 注意这并不是c的源代码。但是这个应该可以说明它是怎么样运行的。当x pi是false 时，逻辑“与”（&）运算后得到的是0，也就是（x-=0）大小没有改变，哈哈完美的等价!我会给读这篇文章的读者留一些关于这个练习。虽然cos 比sin需要多一些运算，但是相比之下貌似也没有更好方法可以让程序更快了。现在我们的最大误差是0.056 ，四项泰勒级数展开式每一次都会有一点点误差。再来看看我们sin函数： 现在是不是不能继续提升精准度了呢？当前的版本已经可以满足大多度sin函数的应用了。但是对一些要求更高一些的程序现在做的还够。仔细观察
8、图像，你会注意到我们的近似值总是比真实值大，当然除了0，pi/2 和 pi。所以我们要做的就是在不改变这些点（0，pi/2 和 pi）的情况下，将函数再“按下去”一些。解决方法是利用抛物线的平方。看起来就像这样： 注意它保持着原来那些关键点，不同的是它比真实的sin函数值更低了。所以我们可以用一个加权的平均值来使两个函数更接近。 Code: Q (4/pi x - 4/pi2 x2) + P (4/pi x - 4/pi2 x2)2 利用Q + P = 1. 你可以灵活的控制绝对误差或相对误差。别急我来告诉你取不同的极限结果时Q，P的值。绝对误差的最佳权值是：Q = 0.775, P = 0.
9、225 ；相对误差的最佳权值是：Q = 0.782，P = 0.218 。让我们来看一下结果的图像。 红线呢？它几乎被绿线完全覆盖了，这足以证明我们的近似十分完美。最大误差是0.001，50倍的提升！这个公式看起来很长，但是括号里面的公式最终得到的值是相同的，也就是说括号里的只需要被计算一次。事实上在原来的基础上只是增加了额外的2次乘法和2次加法就可以得到现在的结果。 先别高兴的太早，我们还要“制造”一个负号出来。我们需要增加一个abs()运算。最终的c代码是： Code: float sine(float x) const float B = 4/pi; const float C = -4
10、/(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项泰勒级数展开式快，更精准！Cos只需要相应的变换一下x就可以了。 (译者注：后面是汇编程序，不翻译了) part2 我选取了最小误差的情况，用as3运行后
11、发现提升了14倍，而且仍然是非常精准。不过你必须直接使用它，不能把它放到一个函数中，因为每调用一次额外的函数调用会削减执行效率，最终你会得到一个比Math.sin() 和 Math.cos()效率更差的结果。 还有这里会用到的三角定理： cos(x) = sin(x + pi/2) cos(x - pi/2) = sin(x) 下载: fastTrig.as. 可以清楚到对比结果，现在你可以用这个替换Math.sin() 和 Math.cos()了 哇哦！几乎是相同的精准度（14倍速度提升） /always wrap input angle to -PI.PIif (x 3.14159265)
12、 x-= 6.28318531;/compute sineif (x 3.14159265) x-= 6.28318531;if (x 0) cos= 1.27323954 * x+ 0.405284735 * x* xelse cos= 1.27323954 * x- 0.405284735 * x* x; High precision sine/cosine (8x faster) /always wrap input angle to -PI.PIif (x 3.14159265) x-= 6.28318531;/compute sineif (x 0) sin= 1.27323954
13、* x+ .405284735 * x* x;if (sin 0) sin= .225 * (sin*-sin- sin) + sin;else sin= .225 * (sin* sin- sin) + sin;else sin= 1.27323954 * x- 0.405284735 * x* x;if (sin 3.14159265) x-= 6.28318531;if (x 0) cos= 1.27323954 * x+ 0.405284735 * x* x;if (cos 0) cos= .225 * (cos*-cos- cos) + cos;else cos= .225 * (c
14、os* cos- cos) + cos;else cos= 1.27323954 * x- 0.405284735 * x* x;if (cos 0) cos= .225 * (cos*-cos- cos) + cos;else cos= .225 * (cos* cos- cos) + cos; Fast and accurate sine/cosine approximation July 18, 2007 on 2:31 pm | In Actionscript | 50 Comments Trigonometric functions are costly operations and
15、 can slow down your application if they are extensively used. There are two reasons why: First, Math.sin() is a function, and thus needs a function call which simple eats up some time. Second, the result is computed with much more precision than you would ever need in most situations. Most often you
16、 just want the periodic wave-like characteristics of the sine or cosine, which can be approximated in various ways. One common way of making it faster is to create a lookup-table by computing the sine at discrete steps and storing the result in an array. For example: var sineTable:Array = ;for (var
17、i:int = 0; i 90; i+) sineTablei = Math.sin(Math.PI/180 * i) Due to the symmetry of the sine wave, its sufficient to compute one quadrant only (0.pi/2), and the other 3/4s of the circle can be computed by shifting and wrapping the input value. The biggest drawback is that the values are stored at a f
18、ixed resolution and so the result is not very accurate. This can be enhanced with linear interpolation: x = 22.5;y = sineTableint(x) + (sineTableint(x + .5) - sineTableint(x) / 2; Much better, but yet the error exists. It also involves accessing array elements which makes the code rather slow. Anoth
19、er technique uses taylor series approximation: sin(x) = x - (x3)/3! + (x5)/5! - (x7)/7! + . Like with the lookup-table, evaluating this term is costly. After searching for alternatives, I finally found a fantastic solution using a quadratic curve which blows everything away in terms of performance a
20、nd accuracy. For a detailed derivation, please follow the link because I wont go into it. I did minor optimizations to figure out what AS3 likes most, and arrived at some code that can be up to 14x faster, while still being very accurate. However, you have to use it directly - do not place the code
21、inside a function, because the additional function call sweeps out the performance gain, and you are left with an approximation that is actually slower compared to a native Math.sin() or Math.cos() call. Also note that cos(x) = sin(x + pi/2) or cos(x - pi/2) = sin(x), so computing the cosine is just
22、 of matter adding pi/2 to the input value. Download source: fastTrig.as. Below is a simple visualization to show you the quality of the approximation. The high precision version can replace the Math.sin() and Math.cos() calls in nearly all situations. Low precision sine/cosine (14x faster) /always w
23、rap input angle to -PI.PIif (x 3.14159265) x -= 6.28318531;/compute sineif (x 3.14159265) x -= 6.28318531;if (x 0) cos = 1.27323954 * x + 0.405284735 * x * xelse cos = 1.27323954 * x - 0.405284735 * x * x; High precision sine/cosine (8x faster) /always wrap input angle to -PI.PIif (x 3.14159265) x -
24、= 6.28318531;/compute sineif (x 0) sin = 1.27323954 * x + .405284735 * x * x;if (sin 0) sin = .225 * (sin *-sin - sin) + sin;else sin = .225 * (sin * sin - sin) + sin;else sin = 1.27323954 * x - 0.405284735 * x * x;if (sin 3.14159265) x -= 6.28318531;if (x 0) cos = 1.27323954 * x + 0.405284735 * x *
25、 x;if (cos 0) cos = .225 * (cos *-cos - cos) + cos;else cos = .225 * (cos * cos - cos) + cos;else cos = 1.27323954 * x - 0.405284735 * x * x;if (cos 0) cos = .225 * (cos *-cos - cos) + cos;else cos = .225 * (cos * cos - cos) + cos; 50 COMMENTS RSS feed for comments on this post. TrackBack URI 1. Nic
26、e find, though I think Id have to be in a real optimization bind before I replaced all my sin/cos calls with all that stuff. :) Comment by Keith Peters July, 18 2007 # 2. Once again, incredible post. There is an error in the link. The good link is :http:/lab.polygonal.de/wp-content/articles/fast_tri
27、g/fastTrig.as :/1.27323954 = 4/pi/0.405284735 =-4/(pi2)/* * low precision sine/cosine */always wrap input angle to -PI.PIif (x 3.14159265) x -= 6.28318531;/compute sineif (x 3.14159265) x -= 6.28318531;if (x 0) cos = 1.27323954 * x + 0.405284735 * x * xelse cos = 1.27323954 * x - 0.405284735 * x * x
28、;/* * high precision sine/cosine */always wrap input angle to -PI.PIif (x 3.14159265) x -= 6.28318531;/compute sineif (x 0) sin = 1.27323954 * x + .405284735 * x * x; if (sin 0) sin = .225 * (sin *-sin - sin) + sin; else sin = .225 * (sin * sin - sin) + sin;else sin = 1.27323954 * x - 0.405284735 *
29、x * x; if (sin 3.14159265) x -= 6.28318531;if (x 0) cos = 1.27323954 * x + 0.405284735 * x * x if (cos 0) cos = .225 * (cos *-cos - cos) + cos; else cos = .225 * (cos * cos - cos) + cos;else cos = 1.27323954 * x - 0.405284735 * x * x; if (cos y) & z; But the code above is a lot slower than the pure
30、if.else statements. Comment by Michael July, 20 2007 # 8. Nice work Michael this one goes a long way back. I think it may have been in Harts Computer Approximations. One thing Ive found helpful when working with compilers that dont support function inlining is to use the C preprocessor or awk or som
31、ething similar to do all my debugging with library calls and then have the tool replace them with inline approximations going into production. Another benefit of having inline code is that you can combine the sin or cos computation with whatever equation its used in and realize additional gains by h
32、aving the trig function result immediately available in a register. regards, - jim Comment by Jim Armstrong July, 20 2007 # 9. hey cool I didnt know there is even a whole book about approximations only ;-) talking about preprocessors, what IDE do you use ? I was trying to integrate the C preprocesso
33、r into FlashDevelop, but failed at configuring make or ant to use it. Comment by Michael July, 21 2007 # 10. Michael: Another classic is Codys A Software Manual for the Elementary Functions. There is a relatively new book by Muller Elementary Functions: Algorithms and Implementation, but I havent re
34、ad it yet. Ive given up on direct integration of other tools to any IDE dont have the bandwidth to work out all the issues, so I run awk or another text processor directly on the .AS files very old school. Best of luck in your continued (and superb) efforts! regards, - jim Comment by Jim Armstrong J
35、uly, 22 2007 # 11. Hey thanks for the reference!Doing University and making a blog isnt easy but shortly some content will popout =D Comment by Eduardo August, 9 2007 # 12. . Fast and accurate sine/cosine approximation More code optimization to do sine/cosine operations without using the Math class.
36、 . Pingback by Blog Archive Links: Mathematical themes with Actionscript polyG August, 26 2007 # 13. . Fast and accurate sine/cosine approximation . Pingback by BLOG Blog Archive as3“ September, 19 2007 # 14. Heres a modified fixed-point version that performs better : http:/blog.haxe.org/entry/26 Co
37、mment by Nicolas November, 12 2007 # 15. Iev tested this in AS2, and its not faster. So this optimisation trick is ONLY for AS3. Just to let you know. Comment by JP January, 8 2008 # 16. Thanks. This post was useful and more accurate than using an MS-Excell trendline function (2nd degree). Made minor changes that should speed things up since addition is cheaper than mult usually. This is the c code:（注：我本人已经亲自测试过OK,性能很好-tang）double sin(double rad)double sine;if (rad 0)sine = rad*(1.27323954 + .405284735 * rad);elsesine = rad*(1.27323954 0.405284735 * rad); if (sine 0)sine = si
