分享

高精度求反正弦、反余弦、反正切、反余切数值算法

 花园图书馆88 2013-09-05

1、反正切算法

在这里以反正切作为文章的开头部分,那是因为反正弦、反余弦、反余切的求值,最终都可以归结为反正切的求解。它们都可以相互推导,在这里鉴于反正切的算法收敛速度快,特以此作为最基本的函数。

这里我们记x的反正切形式为Arctan(x),由于Arctan(-x)=-Arctan(x),所以我们这里考虑的x是大于0的情形,小于0的情形应该很容易处理吧,这里不再多说。根据泰勒展开有

Arctan(x)=x-x^3/3+x^5/5+……+(-1)^(n+1)*x^(2n+1)/(2n+1)+……

=x*(1-x^2/3+x^4/5+……+(-1)^(n+1)*x^(2n)/(2n+1)+……)

如果我们设y(n)=1-x^2/3+x^4/5+……+(-1)^(n+1)*x^(2n)/(2n+1)+……

则Arctan(x)=x*y(n)

在这里我们对y(n)可以很容易构造递推式子

y(n)=1/(2n-1)-n^2*y(n+1)

当n越大,x*y(n)就越接近我们需要的结果,这样就有如下的算法

【参考《实用数值计算方法》.甄西丰.清华大学出版社.2006.76页】

算法1.1

(1)、设y(n+1)=1/(2n+1)

(2)、对于k=n,n-1,n-2,……,1进行计算y(k)=1/(2k-1)-x^2*y(k+1)

(3)、最后返回x*y(1)的值

其中误差=x^2n/(2n+1)

优化算法:

从上面的误差可以看出,当|x|<1时才收敛稍微快一点,那有没有办法让x尽可能的小呢,这里是有办法的.看下面一个公式

Arctan(x)=Arctan(y)+Arctan((x-y)/(1+x*y))

即我们计算Arctan(x)只要我们计算Arctan(y)与Arctan((x-y)/(1+x*y))的值就可以了.利用这个公式我们可以构造一个很酷的算法.

Arctan(x)

{

if (x<0)

{ 返回-Arctan(-x)}

elseIf(x<0.25)

{ 直接使用上面的算法1.1进行求解Arctan(x) ;}

elseIf(x<0.75)

{ 返回Arctan(x)=Arctan(0.5)+Arctan((x-0.5)/(1+x*0.5));}

elseIf(x<1)

{ 返回Arctan(x)=Arctan(0.75)+Arctan((x-0.75)/(1+x*0.75));}

elseIf(x<2)

{ 返回Arctan(x)=Arctan(1.5)+Arctan((x-1.5)/(1+x*1.5));}

else

{ 返回Arctan(x)=Arctan(4)+Arctan((x-4)/(1+x*4));;}

}

上面的Arctan(0.75)、Arctan(1.5)、Arctan(4)的值可以提前先计算出来,上面判断分别用了0.25、0.75、1、2等值,当然你可以使用其它的值.在这里,使用上面的Arctan(x)算法,计算中最多只调用算法1.1一次,就可以得到想要的结果。根据上面这个,在算法1.1中,当n=100的时候,最差的情况=0.25^200/(200+1)≈1.92666264420364E-123即可保证到小数点后120多位,当然实际应该更多,如果你n取更大,判断的值比0.25更小,则计算的有效位数也就更多.这里不再多说.

2、反余切算法

我们知道Tan(x)=a,Cot(x)=1/a,那么Arccot(a)=Arctan(1/a),所以求Arccot(x)的时候直接调用上面的Arctan(x)算法即可.

3、反正弦算法

我们知道在复数领域内参看《复数的基本运算》

Arcsin(x)=-i*Ln[x*i±(1-x^2)^0.5]=-i*[ln|x*i±(1-x^2)^0.5|+Arg(x*i±(1-x^2)^0.5)*i]=Arg(x*i±(1-x^2)^0.5)

在这里我们就可以构造如下的算法

算法3.1

Arcsin(x)

{

If (x<0)

{ return -Arcsin(-x);}

ElseIf(x<1)

{return Arctan(x/(1-x^2)^0.5);}

Else

{ Return π/2;}

}

4、反余弦算法

我们知道在复数领域内参看《复数的基本运算》

Arccos(x)=-i*Ln[x±(x^2-1)^0.5]=-i*[ln|x±(x^2-1)^0.5|+Arg(x±(x^2-1)^0.5)*i]=Arg(x±(x^2-1)^0.5)=Arg(x±(1-x^2)^0.5*i)

在这里我们就可以构造如下的算法

算法4.1

Arccos(x)

{

If (x=-1)

{Return π;}

ElseIf(x<0)

{ Return π-Arctan((1-x^2)^0.5/(-x));}

ElseIf(x=0)

{return π/2;}

Elseif (x<1)

{ Return Arctan((1-x^2)^0.5/x;}

Else

{ Return 0;}

}

注意:

(1)上面求反正弦与反余弦的时候都涉及到求解Arctan(x^0.5)这样的问题,我们知道对一个大数x进行开方的成本是很高的,至少目前我接触到的算法效率都很低.所以我们能否找到一种关系,比如我现在已经计算出Arctan(x),Arccot(x),Sin(x),Cos(x),Tan(x)的值,我能否通过这些计算出来的值来推导Arctan(x^0.5)的值呢?拭目以待吧.如果您有幸看到这里,并且你如果能推导出具体的公式,请给我留言一份,不甚感激!

(2)当然,在上面求解反正弦与反余弦的时候可以参考【《实用数值计算方法》.甄西丰.清华大学出版社.2006.70-72页】提到的泰勒展开,但是文章中我为什么没有用他的算法,因为当x趋近0的时候收敛才会快,趋近1的时候就相当慢了.当然如果能找到类似Arctan(x)=Arctan(y)+Arctan((x-y)/(1+x*y))这样的公式,那就不用反正切来计算反余切了.这样如果能用更小的数来计算更大的数,那就是相当漂亮的.但是目前能力有限,还没有推导出来,有空专门去推导一下.

源代码:反正切反余切,反正弦与反余弦

    本站是提供个人知识管理的网络存储空间,所有内容均由用户发布,不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息,谨防诈骗。如发现有害或侵权内容,请点击一键举报。
    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多