高精度反正切函数的实现1

本帖最后由 落叶 于 2016-6-5 09:16 编辑 经过了几天艰难的查资料，对高精反正切函数实现并运行成功，分亨一下， 主要资料参考：http://blog.163.com/shikang999@1 ... 962012426103454943/ 因为加入了我的改进，和我自已的方法，所以我在原创中发布： 这里反正切函数的根本实现方法还是泰勒展开式，但它不能单独采用，必须使自变量X趋进这个函数曲线的平滑度最高的点，这个点的收敛速度最快， 反正切函数arctan(x)的这个点是0，反正弦函数arcsin(x),反余弦函数arccos(x)的最快收敛点也是零， 我看了N遍的反正弦三角函数方面的公式，并没有找到使X缩小的公式，不过世慷（上述参考资料的发布者）找到了一个方法： 通过这几个公式实现反正切函数: 根本公式为反正切公式：Arctan(x)=Arctan(y)+Arctan((x-y)/(1+x*y)) 其它公式，Arccot(a)=Arctan(1/a)，arcsin(x)=Arctan(x/(1-x^2)^0.5)，arccos(x)=π-Arctan((1-x^2)^0.5/(-x))    （x<0）， arccos(x)=Arctan((1-x^2)^0.5/x     （0<x<1） 世慷的具体方法如下：为了方便讨论，把它的源语句复制过来： 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(4)，Arctan(1.5)，Arctan(0.75)，Arctan(0.5) 但它只能对X进行一次优化，效率较低，另外自变量优化后范围较分散，不利于后面的泰勒运算统一运算级数，会有很大的冗余运算，再就是四个常数取的太随意。 我的方法如下： 取Arctan(1)，Arctan(0.1)，Arctan(0.01)，Arctan(0.001)，Arctan(0.0001)        预算出来 第一个取值1很重要，因为所有大于1的数经过公式分解后都小于1，                 temp = string1                    ‘string1为自变量X                 j = 1                                    ’g01是一个高精数组，g01(1)=1，g01(2)=0.1......g01(5)=0.0001                 Do                     Do                         x = myCompare(temp, g01(j))               '大数比较函数，大于返回2,小于返回1，等于返回0                         If x < 2 Then                             Exit Do                         End If                             temp1 = myMINUS(temp, g01(j))                      ‘大数减，vb6.0不支持运算符重载，只能这样写了                             temp2 = myMULTIPLY(temp, g01(j))                 ’大数乘                             temp2 = myADD(g1, temp2)                             ‘大数加                             temp = myEXCEPT(temp1, temp2)                    ’大数除                              tempCs = myADD(tempCs, fzqb(j))                    ‘这步是关键，例：假如X经过了1分解1 次，0.1分解5次.......                     Loop                                                                           'tempCs的值就是Arctan(1)*1+Arctan(0.1)*5.......                     j = j + 1                 Loop Until j = 6                 temp = myatTan(temp)                                                  '对分解后的x进行泰勒运算                 temp = myADD(temp, tempCs) 上述方法可以不断使X变小，使X小于0.0001，对大于1的值，几个常数的利用率接近百分之百。 现在还有一个问题：开始我的常数是这样的，1，    0.001 ，     0.00001 ，0.0000001, 0.000000001  X分解的速度远慢于泰勒级数的运算速度，所以我取了1,0.1,0.001....，但现在X分解的速度又快于泰勒级数的运算速度，但我已经没有精力再测试了，两者速度大致相等效率最高。 有兴趣的朋友可以测测各个常数之间的数距多少合适？ 因为上面的X优化，现泰勒运算的级数采用（精度除以4），也可获得正确结果，世慷的算法需要运算的级数（精度乘以2-3） 泰勒计算公式为：Arctan(x)=x-x^3/3+x^5/5+……+(-1)^(n+1)*x^(2n+1)/(2n+1) 修改为：Arctan(x)=x-(x*x^2（Gb/3-x^2（Gb/5-x^2(Gb/7-……)))/Gb  (-1)^(n+1)*x^(2n+1)/(2n+1) Gb为1  3   5   7  .......的公倍数 经过这个公式优化后，万位精度运算大概获得4-5秒的加速 现在我的程序的万位精度运算耗时26秒左右， 但BTCAL疯狂计算器只需要2-3秒，太打击我了，也不知这个函数是他自已写的，还是调用的大数库？不过他也有BUG，在万位Arctan（1.1）运算时，用时高达6分50秒，接近1的值都出现运算困难问题。 我的程序所有值都慢，都是一个速度，26秒左右。 我就说到这里，欢迎网友们批评，指正。 现这个方法我已没有采用了，得到BTCAL疯狂计算器作者的帮助，采用了和他一样的算法，并用我推算的公式优化了接近1的值运算困难问题。