找回密码
 欢迎注册
查看: 14206|回复: 4

[原创] 高精度反三角函数的实现

[复制链接]
发表于 2016-6-5 09:13:37 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
经过了几天艰难的查资料,对高精反三角函数实现并运行成功,分亨一下,
主要资料参考: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秒左右。
我就说到这里,欢迎网友们批评,指正。

补充内容 (2016-6-7 11:41):
前面所说(泰勒运算的级数)并不是真实的级数,为了编程方便,我把运算最后一项的除数作为级数,真实运算级数是我上述说法的一半左右。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-6-8 07:22:20 | 显示全部楼层
反三角函数可以以三角函数为基础,使用牛顿迭代法求解。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-6-8 14:16:11 | 显示全部楼层
只是不知牛顿迭代式是怎么样的?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-6-8 16:39:50 | 显示全部楼层
如果知道反三角函数的导数,牛顿迭代公式不难求出来。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-6-8 16:57:29 | 显示全部楼层
用 算术-几何平均求复数对数吧
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-11-21 23:50 , Processed in 0.024611 second(s), 17 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表