三角函数逼近快速算法正余弦.docx
- 文档编号:1438746
- 上传时间:2023-05-01
- 格式:DOCX
- 页数:23
- 大小:70.38KB
三角函数逼近快速算法正余弦.docx
《三角函数逼近快速算法正余弦.docx》由会员分享,可在线阅读,更多相关《三角函数逼近快速算法正余弦.docx(23页珍藏版)》请在冰点文库上搜索。
三角函数逼近快速算法正余弦
三角函数逼近快速算法(正余弦)
原文出自:
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函数。
有时候我们并需要过高的精度,这时C语言中自带的三角函数(sinf()和cosf()f)计算的精度超出了我们所需要的精度要求,所以其效率很低。
我们真正需要的是间于精度和效率的一个折中的方案。
众所周知的取近似值的方法是:
泰勒级数(和著名的马克劳林级数)
代码是:
x-1/6x^3+1/120x^5-1/5040x^7+...
我们绘制了下图:
绿线是标准的sin函数,红线是4项泰勒级数展开式。
这个近似值的效果看起来还不错,但是如果你仔细观察后会发现
它在pi/2之前的效果还是很好的,但是超过了pi/2后就开始快速偏离标准sin。
它在pi处的值比原来的0多了0.075.用这个方法来模拟波动忽动忽停,看起来很呆板,这个效果当然不是我们想要的。
我们可以继续增加泰勒级数项的个数来减小误差,但是这将导致我们的公式非常的冗长。
用4项的泰勒级数展开式需要我们进行7次乘法和3次加法来完成。
泰勒级数不能同时满足我对精度和效率的要求。
刚刚近似如果能满足sin(pi)=0就好了。
从上图我还可以发现另一件事:
这个曲线看起来很象抛物线!
所以我们来寻找一个尽可能和sin接近的抛物线(公式)。
抛物线的范式方程是:
A+Bx+Cx^2.这个公式可以让们控制三个自由度。
显然我们需要其满足sine(0)=0,sine(pi/2)=1andsine(pi)=0.这样我们就得到了3个等式。
A+B0+C0^2=0
A+Bpi/2+C(pi/2)^2=1
A+Bpi+Cpi^2=0
解得:
A=0,B=4/pi,C=-4/pi^2.我们的抛物线诞生啦!
貌似这个的误差看起来比泰勒级数还要遭。
其实不是的!
这种方法的最大误差是0.056.(译者:
而且这种近似值没有误差积累)而且这个近似值的绘制出的波动是光滑的,而且只需要3次乘法和一次加法。
不过它还不够完美。
下图是[-pi,pi]之间的图像:
显然我们至少需要它在1个完整的周期内都符合我们要求。
但是我们可以看出,我们缺少的另一半是原抛物线的一个映射。
它的公式是:
4/pix+4/pi^2x^2。
所以我们可以直接这样写:
Code:
if(x>0){y=4/pix-4/pi^2x^2;}else{y=4/pix+4/pi^2x^2;}
添加一个条件分支不是一个好的方法。
它会让程序渐渐的变慢。
但是观察一下我们模拟的和标准的图像是多么的接近啊!
观察上面两式子,只是中间的一项正负号不同,我的第一个单纯的想法是可以提取x的正负号来消除分支,即使用:
x/abs(x)。
除法的代价是非常大的,我们来观察一下现在的公式:
4/pix-x/abs(x)4/pi^2x^2。
将除法化简后我们得到:
4/pix-4/pi^2xabs(x).所以只需要多一步运算我们就得到了我们需要的sin逼近值!
下图是结果
接下来我们要考虑cos。
有基础的三角公理可以知道:
cos(x)=sin(pi/2+x).把x多加一个pi/2就可以搞定了?
事实上它的某一部分不是我们期望得到的。
我们需要做的就是当x>pi/2时“跳回”。
这个可以由减去2pi来实现。
Code:
x+=pi/2;
if(x>pi)//Originalx>pi/2{x-=2*pi;//Wrap:
cos(x)=cos(x-2pi)}
y=sine(x);
又出现了一个分支,我们可以用逻辑“与”来消除它,像是这样:
x-=(x>pi)&(2*pi);
Code:
x-=(x>pi)&(2*pi);
注意这并不是c的源代码。
但是这个应该可以说明它是怎么样运行的。
当x>pi是false时,逻辑“与”(&)运算后得到的是0,也就是(x-=0)大小没有改变,哈哈完美的等价!
我会给读这篇文章的读者留一些关于这个练习。
虽然cos比sin需要多一些运算,但是相比之下貌似也没有更好方法可以让程序更快了。
现在我们的最大误差是0.056,四项泰勒级数展开式每一次都会有一点点误差。
再来看看我们sin函数:
现在是不是不能继续提升精准度了呢?
当前的版本已经可以满足大多度sin函数的应用了。
但是对一些要求更高一些的程序现在做的还够。
仔细观察图像,你会注意到我们的近似值总是比真实值大,当然除了0,pi/2和pi。
所以我们要做的就是在不改变这些点(0,pi/2和pi)的情况下,将函数再“按下去”一些。
解决方法是利用抛物线的平方。
看起来就像这样:
注意它保持着原来那些关键点,不同的是它比真实的sin函数值更低了。
所以我们可以用一个加权的平均值来使两个函数更接近。
Code:
Q(4/pix-4/pi^2x^2)+P(4/pix-4/pi^2x^2)^2
利用Q+P=1.你可以灵活的控制绝对误差或相对误差。
别急我来告诉你取不同的极限结果时Q,P的值。
绝对误差的最佳权值是:
Q=0.775,P=0.225;相对误差的最佳权值是:
Q=0.782,P=0.218。
让我们来看一下结果的图像。
红线呢?
它几乎被绿线完全覆盖了,这足以证明我们的近似十分完美。
最大误差是0.001,50倍的提升!
这个公式看起来很长,但是括号里面的公式最终得到的值是相同的,也就是说括号里的只需要被计算一次。
事实上在原来的基础上只是增加了额外的2次乘法和2次加法就可以得到现在的结果。
先别高兴的太早,我们还要“制造”一个负号出来。
我们需要增加一个abs()运算。
最终的c代码是:
Code:
floatsine(floatx)
{
constfloatB=4/pi;
constfloatC=-4/(pi*pi);
floaty=B*x+C*x*abs(x);
#ifdefEXTRA_PRECISION//constfloatQ=0.775;
constfloatP=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倍速度提升)
//alwayswrapinputangleto-PI..PI
if(x<-3.14159265)
x+=6.28318531;
else
if(x>3.14159265)
x-=6.28318531;
//computesine
if(x<0)
sin=1.27323954*x+.405284735*x*x;
else
sin=1.27323954*x-0.405284735*x*x;
//computecosine:
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;
}
Highprecisionsine/cosine(~8xfaster)
//alwayswrapinputangleto-PI..PI
if(x<-3.14159265)
x+=6.28318531;
else
if(x>3.14159265)
x-=6.28318531;
//computesine
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;
}
//computecosine:
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;
}
Fastandaccuratesine/cosineapproximation
July18,2007on2:
31pm|InActionscript|50Comments
Trigonometricfunctionsarecostlyoperationsandcanslowdownyourapplicationiftheyareextensivelyused.Therearetworeasonswhy:
First,Math.sin()isafunction,andthusneedsafunctioncallwhichsimpleeatsupsometime.Second,theresultiscomputedwithmuchmoreprecisionthanyouwouldeverneedinmostsituations.
Mostoftenyoujustwanttheperiodicwave-likecharacteristicsofthesineorcosine,whichcanbeapproximatedinvariousways.Onecommonwayofmakingitfasteristocreatealookup-tablebycomputingthesineatdiscretestepsandstoringtheresultinanarray.Forexample:
varsineTable:
Array=[];
for(vari:
int=0;i<90;i++)
{
sineTable[i]=Math.sin(Math.PI/180*i)
}
Duetothesymmetryofthesinewave,it'ssufficienttocomputeonequadrantonly(0..pi/2),andtheother3/4'softhecirclecanbecomputedbyshiftingandwrappingtheinputvalue.Thebiggestdrawbackisthatthevaluesarestoredatafixedresolutionandsotheresultisnotveryaccurate.Thiscanbeenhancedwithlinearinterpolation:
x=22.5;
y=sineTable[int(x)]+(sineTable[int(x+.5)]-sineTable[int(x)])/2;
Muchbetter,butyettheerrorexists.Italsoinvolvesaccessingarrayelementswhichmakesthecoderatherslow.Anothertechniqueusestaylorseriesapproximation:
sin(x)=x-(x^3)/3!
+(x^5)/5!
-(x^7)/7!
+...
Likewiththelookup-table,evaluatingthistermiscostly.
Aftersearchingforalternatives,Ifinallyfoundafantasticsolutionusingaquadraticcurvewhichblowseverythingawayintermsofperformanceandaccuracy.Foradetailedderivation,pleasefollowthelinkbecauseIwon'tgointoit.
IdidminoroptimizationstofigureoutwhatAS3likesmost,andarrivedatsomecodethatcanbeupto14xfaster,whilestillbeingveryaccurate.However,youhavetouseitdirectly-donotplacethecodeinsideafunction,becausetheadditionalfunctioncallsweepsouttheperformancegain,andyouareleftwithanapproximationthatisactuallyslowercomparedtoanativeMath.sin()orMath.cos()call.Alsonotethatcos(x)=sin(x+pi/2)orcos(x-pi/2)=sin(x),socomputingthecosineisjustofmatteraddingpi/2totheinputvalue.
Downloadsource:
fastTrig.as.
Belowisasimplevisualizationtoshowyouthequalityoftheapproximation.ThehighprecisionversioncanreplacetheMath.sin()andMath.cos()callsinnearlyallsituations.
Lowprecisionsine/cosine(~14xfaster)
//alwayswrapinputangleto-PI..PI
if(x<-3.14159265)
x+=6.28318531;
else
if(x>3.14159265)
x-=6.28318531;
//computesine
if(x<0)
sin=1.27323954*x+.405284735*x*x;
else
sin=1.27323954*x-0.405284735*x*x;
//computecosine:
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;
}
Highprecisionsine/cosine(~8xfaster)
//alwayswrapinputangleto-PI..PI
if(x<-3.14159265)
x+=6.28318531;
else
if(x>3.14159265)
x-=6.28318531;
//computesine
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;
}
//computecosine:
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;
}
50COMMENTS»
RSSfeedforcommentsonthispost.TrackBackURI
1.Nicefind,thoughIthinkI’dhavetobeinarealoptimizationbindbeforeIreplacedallmysin/coscallswithallthatstuff.:
)
CommentbyKeithPeters—July,182007#
2.Onceagain,incrediblepost.
Thereisanerrorinthelink.Thegoodlinkis:
http:
//lab.polygonal.de/wp-content/articles/fast_trig/fastTrig.as:
//1.27323954 =4/pi
//0.405284735=-4/(pi^2)
/*********************************************************
*lowprecisionsine/cosine
*********************************************************/
//alwayswrapinputangleto-PI..PI
if(x<-3.14159265)
x+=6.28318531;
else
if(x> 3.14159265)
x-=6.28318531;
//computesine
if(x<0)
sin=1.27323954*x+.405284735*x*x;
else
sin=1.27323954*x-0.405284735*x*x;
//computecosine:
sin(x+PI/2)=cos(x)
x+=1.57079632;
if(x> 3.14159
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三角函数 逼近 快速 算法 余弦