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))这样的公式,那就不用反正切来计算反余切了.这样如果能用更小的数来计算更大的数,那就是相当漂亮的.但是目前能力有限,还没有推导出来,有空专门去推导一下. 源代码:反正切反余切,反正弦与反余弦 |
|