高精度计算反正弦反余弦反正切反余切算法

shikang999

[这个贴子最后由shikang999在 2012/06/10 03:29pm 第 1 次编辑] 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;} } 问题 上面的反正弦与反余弦的算法应该还有优化的余地,我想问下陆教授,能否像公式Arctan(x)=Arctan(y)+Arctan((x-y)/(1+x*y))一样找到一个Arcsin(x)=Arcsin(y)+Arcsin(z)或者Arcsin(x)=Arccos(y)+Arcsin(z)或者Arcsin(x)=Arcsin(y)+Arctan(z)这样的公式呢?如果能找到的话,也许这样计算量会减小很多。 例子,下面是我编程计算的结果[\b] Arctan(1.23456789123456789123456789) = "0.8899874939173517418007233056028177395094603517310346085764627411189975002207321929061547891934656229545228592405633992500819082078686173218162428898059394528442704538294231492204183061206689827196142356149669443370434987090492989000212486347187469000368907156877608031324819612128086872467478982041579370567830903027758810382550171180053155721534269222210921218128923103591685379139588787488837804163235657630499639428197045431009114098768853603455282160499155147715360419989739777682162515"