posts - 72,  comments - 34,  trackbacks - 0

 看过第一篇的文章后,大呼过瘾!原文作者的思路非常简捷,有趣,偶英语比较差,欢迎指正,废话不多说看文章

原文出处:

http://www.devmaster.net/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 > pifalse 时,逻辑“与”(&)运算后得到的是0,也就是(x-=0)大小没有改变,哈哈完美的等价!我会给读这篇文章的读者留一些关于这个练习。虽然cos 比sin需要多一些运算,但是相比之下貌似也没有更好方法可以让程序更快了。现在我们的最大误差是0.056 ,四项泰勒级数展开式每一次都会有一点点误差。再来看看我们sin函数:

 

现在是不是不能继续提升精准度了呢?当前的版本已经可以满足大多度sin函数的应用了。但是对一些要求更高一些的程序现在做的还够。

仔细观察图像,你会注意到我们的近似值总是比真实值大,当然除了0,pi/2  pi。所以我们要做的就是在不改变这些点(0,pi/2  pi)的情况下,将函数再“按下去”一些。解决方法是利用抛物线的平方。看起来就像这样:

注意它保持着原来那些关键点,不同的是它比真实的sin函数值更低了。所以我们可以用一个加权的平均值来使两个函数更接近。

Code:

Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2

利用Q + P = 1. 你可以灵活的控制绝对误差或相对误差。别急我来告诉你取不同的极限结果时Q,P的值。绝对误差的最佳权值是:Q = 0.775, P = 0.225 ;相对误差的最佳权值是:Q = 0.782P = 0.218 。让我们来看一下结果的图像。 

 

红线呢?它几乎被绿线完全覆盖了,这足以证明我们的近似十分完美。最大误差是0.001,50倍的提升!这个公式看起来很长,但是括号里面的公式最终得到的值是相同的,也就是说括号里的只需要被计算一次。事实上在原来的基础上只是增加了额外的2次乘法和2次加法就可以得到现在的结果。

先别高兴的太早,我们还要“制造”一个负号出来。我们需要增加一个abs()运算。最终的c代码是: 

Code:

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项泰勒级数展开式快,更精准!Cos只需要相应的变换一下x就可以了。

(译者注:后面是汇编程序,不翻译了)

part2

 我选取了最小误差的情况,用as3运行后发现提升了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..PI
if (x< -3.14159265)
    x+= 6.28318531;
else
if (x>  3.14159265)
    x-= 6.28318531;

//compute sine
if (x< 0)
    sin= 1.27323954 * x+ .405284735 * x* x;
else
    sin= 1.27323954 * x- 0.405284735 * x* x;

//compute cosine: sin(x + PI/2) = cos(x)
x+= 1.57079632;
if (x>  3.14159265)
    x-= 6.28318531;

if (x< 0)
    cos= 1.27323954 * x+ 0.405284735 * x* x
else
    cos= 1.27323954 * x- 0.405284735 * x* x;
}

High precision sine/cosine (~8x faster)

//always wrap input angle to -PI..PI
if (x< -3.14159265)
    x+= 6.28318531;
else
if (x>  3.14159265)
    x-= 6.28318531;

//compute sine
if (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< 0)
        sin= .225 * (sin*-sin- sin) + sin;
   else
        sin= .225 * (sin* sin- sin) + sin;
}

//compute cosine: sin(x + PI/2) = cos(x)
x+= 1.57079632;
if (x>  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< 0)
        cos= .225 * (cos*-cos- cos) + cos;
   else
        cos= .225 * (cos* cos- cos) + cos;
}

0
0
(请您对文章做出评价)
« 上一篇:【转】【基础】AS3.0里的基元数据类型
» 下一篇:
posted on 2009-03-20 19:41 sungo 阅读(3382) 评论(19)  编辑 收藏 网摘

FeedBack:
2009-03-20 20:23 | 无聊之人[未注册用户]
在游戏开发中这种公式应该很多,不需要太精确的数学公式,为了效率大量简化,极端时用sin表(在没有计算器时,就是靠查表的)来直接查表
  回复  引用    
#2楼[楼主]
2009-03-20 20:44 | 灵药      
--引用--------------------------------------------------
无聊之人: 在游戏开发中这种公式应该很多,不需要太精确的数学公式,为了效率大量简化,极端时用sin表(在没有计算器时,就是靠查表的)来直接查表
--------------------------------------------------------
受教了,我是第一次在flash里面看到这个“新奇”东西,所以很是兴奋;而且作者的原文真的是挺不错,我文笔不好翻不出味道来

  回复  引用  查看    
2009-03-20 23:12 | 徐少侠      
这个不顶更待何时

貌似我的数学早还给大学老师了

所剩无几了

惭愧啊

  回复  引用  查看    
#4楼[楼主]
2009-03-20 23:18 | 灵药      
@徐少侠
-引用--------------------------------------------------
徐少侠: 这个不顶更待何时

貌似我的数学早还给大学老师了

所剩无几了

惭愧啊
--------------------------------------------------------
与君共勉

  回复  引用  查看    
2009-03-21 00:08 | cutepig      
不错不错,不过我觉得在0~pi/2区间用泰勒展开,然后将需要计算的角度换算到该区间在计算应该精度比抛物线高~
  回复  引用  查看    
#6楼[楼主]
2009-03-21 05:54 | 灵药      
@cutepig
恩,不过会额外的添加几个程序分支,不太容易全部都消除了,这样可能会降低效率。但是你说的也很有意思,我会试试的,谢谢你的留言。

  回复  引用  查看    
2009-03-21 13:21 | winter-cn[未注册用户]
分段做二次插值就好了 没必要绕这么大圈子
  回复  引用    
2009-03-21 13:31 | winter-cn[未注册用户]
泰勒级数公式拿到以后也可以优化的
这个有些做法很巧妙 但是缺乏合理性 最终搞出来的还是一个4次函数
还不如把4阶泰勒级数优化一下计算方法

  回复  引用    
#9楼[楼主]
2009-03-21 15:40 | 灵药      
@winter-cn
--引用--------------------------------------------------
winter-cn: 泰勒级数公式拿到以后也可以优化的
这个有些做法很巧妙 但是缺乏合理性 最终搞出来的还是一个4次函数
还不如把4阶泰勒级数优化一下计算方法
--------------------------------------------------------
谢谢你的留言,不过我不太懂你所指的“合理性”是哪一方面的合理性?
你提出的观点我记下了,过一段可能会再出一篇相关的文儿。

  回复  引用  查看    
2009-03-21 15:51 | lookof      
略略看过一点,果然名不虚传。这个想法的确很巧妙的,让偶觉得数学即使科学严肃的,又同时是可以“耍花招”的,很有意思。
等我研究优化时再好好看看。顶~

  回复  引用  查看    
2009-03-21 18:41 | 风过 无痕      
更像是雕虫小技,很多这样的函数拟合,就是多个函数叠加,如果一个函数误差较大,不能满足要求,再取另外一个符号相反的函数来抵消。不仅三角函数,各种没确定解析式的函数拟合...多试几次,加上经验,没什么高深的
  回复  引用  查看    
2009-03-21 18:46 | winter-cn[未注册用户]
@灵药
从二次函数到这个
Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2
我没有看到合理的理由和清晰的思路

从最后这个式子的形状看 这是一个带参数的4次函数的特例 为什么它能接近三角函数 其它4次函数能不能更接近或者运算量更小? 这个地方似乎完全是凭主观感觉写的 所以我说思路不自然 缺乏合理性

  回复  引用    
#13楼[楼主]
2009-03-22 05:40 | 灵药      
@lookof
哈 虽然不严谨,但是我第一次觉得数学这样有意思

@风过 无痕
这些东西可能是您所说的确是存在于“投机取巧”+“经验”之上,但是我第一次看这篇文章感觉是:一种莫名的兴奋,感觉特别有意思,数学也可以这样(虽然不严谨,是大伙提醒我才知道的),菜鸟日志,见笑了。
嘿嘿 另外能否解答一下这个问题呢?我实在是想不出来,数学不太好^_^
--引用--------------------------------------------------
winter-cn: @灵药
从二次函数到这个
Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2
我没有看到合理的理由和清晰的思路

从最后这个式子的形状看 这是一个带参数的4次函数的特例 为什么它能接近三角函数 其它4次函数能不能更接近或者运算量更小? 这个地方似乎完全是凭主观感觉写的 所以我说思路不自然 缺乏合理性
--------------------------------------------------------

  回复  引用  查看    
#14楼[楼主]
2009-03-22 05:49 | 灵药      
@winter-cn
--引用--------------------------------------------------
winter-cn: @灵药
从二次函数到这个
Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2
我没有看到合理的理由和清晰的思路

从最后这个式子的形状看 这是一个带参数的4次函数的特例 为什么它能接近三角函数 其它4次函数能不能更接近或者运算量更小? 这个地方似乎完全是凭主观感觉写的 所以我说思路不自然 缺乏合理性
--------------------------------------------------------
嘿嘿 最初接触这篇文章时,光顾着傻乐文章的趣味性了,当时看的时候下意识的感觉在哪见过所以就默认对了。这个数学问题说实话我暂时还真不知道,而且也没发现 ⊙﹏⊙b汗 我会去试试找出答案的,再次感谢您的留言。

  回复  引用  查看    
2009-03-22 10:41 | 蛙蛙池塘      
这个帖子脑袋竟然不过来凑热闹,呵呵。
  回复  引用  查看    
2009-03-23 12:22 | saimk2[未注册用户]
呃……我觉得还是查表吧
虽然这方法不是特别实用,但是很帅,真的。

  回复  引用    
#17楼[楼主]
2009-03-23 14:27 | 灵药      
@saimk2
哈哈 说的是,不过这种方法很有趣啊^_^

  回复  引用  查看    
2009-05-03 23:26 | jub[未注册用户]
你好,是否可以直接分享完整C程式
謝謝

  回复  引用    
#19楼[楼主]
2009-05-04 16:51 | sungo      
@jub
你好,其实我翻译这篇文章是因为作者的方法和思路让我觉得很兴奋很有趣,c的代码我没有写过,也仅仅是理解,我主要关注的是as3的实现。看过教程我相信写出应该不难,可能是我的翻译有些晦涩,有兴趣的可以去看看原文:http://www.devmaster.net/forums/showthread.php?t=5784

  回复  引用  查看    
<2009年3月>
22232425262728
1234567
891011121314
15161718192021
22232425262728
2930311234

搜索

 

常用链接

我的标签

随笔档案

link

最新评论

阅读排行榜

评论排行榜