找回密码
 欢迎注册
楼主: 落叶

[原创] 基于泰勒展开式的高精三角函数实现

[复制链接]
发表于 2016-4-14 18:27:39 | 显示全部楼层
liangbch 发表于 2016-4-14 18:23
一楼的算法就是 秦九韶算法。

差别是除以阶乘,是每次小乘法最后一次大除法还是每次小除法。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-16 10:20:49 | 显示全部楼层
本帖最后由 落叶 于 2016-4-16 16:58 编辑
zeroieme 发表于 2016-4-14 17:24
是不是小除法比一次大除法慢?

x*(1-(x^2/6)*(1-(x^2/20)*(1-(x^2/42)*(1-x^2/72*(……)))))


这个式子还真是巧妙,竟然把里面的大除法都换成了小除法,肯定会有很大的提速,我的第一个方法只是规避了x的重复运算,和阶乘的重复运算,你这个式子具有同样的功能,并改善了除法,效率更高.
多次小除法和一次大除法,我认为运算次数较多时小除法要慢,小除法并不能充分利用硬件性能,而大除法除了更充分利用硬件外,还有如FFT之类神器加速.
sinx= x- x*x^2(5*7*9 - x^2(7*9 - x^2(9- x^2)))/9!!,这个式子我简单算了一下,结果有点不一致.
我的第二个方法我觉得并不只是多次小除法和一次大除法的问题,当运算项数较多时,可以近似认为消除了除法,所以才有了(针对第一个方法)十倍的提升,所以也可是说这个提速和秦九韶算法关系不大,在把分子多项式化简时,我的思路是这样想的,因为每项之间有一个x^2的增量,我想把它分解出来,看看能不能简化编程,结果最后和秦九韶算法效果一样,这只是一个巧合,不过秦九韶算法真的不错!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-16 16:34:48 | 显示全部楼层
本帖最后由 落叶 于 2016-4-16 18:47 编辑

经过几天的学习,查资料,对sin()的高精算法又进一步改进,分享一下,欢迎吧友们指正。
还是针对sin()函数的泰勒展开式进行计算:
到改进公式:sinx=x−x⋅x2(4⋅5⋅6⋅7⋅8⋅9−x2(6⋅7⋅8⋅9−x2(8⋅9−x2)))/9!
倍角公式:       sin3x=3sinx-4(sinx)^3
倍角公式的算法不是我自已想出来的,是一篇论文上的,他简单分析了一下用倍角公式分解角x,下面的实现是我自已分析的。
倍角公式理论上可以无限把角x分小,而无限趋近于零,在设定精度减少泰勒级数,但这里没有两全其美的事,随着倍角公式的细分,细分级数增加,泰勒公式计算结果返回后的运算量也在急剧增加,我粗略调试了一下,
针对一万精度,倍角公式细分60次,泰勒级数设为精度除以30,运算速度相对我的算法二,速度提高8倍,现在的速度是一万精度sinx,用时13秒左右,(因为里面的关键乘法是我自已写的,其实乘法算法已经采有了比较快的算法,但由于各种原因,速度比BTCAL疯狂计算器v2.5慢10倍,以后有空了再改调用GMP的乘法),下面是流程:
                x = x MOD  2pi                           '弧度化简,x \2pi,取余数,余数为角x,  sinx值不变,x 现在小于2pi
                x = x MOD  pi                            '弧度化简,x \pi,取余数,余数为角x,  sinx绝对值不变, x 现在小于pi
                使  x 小于pi/2                            ‘这一步有点不好表达,咱用中文编程                 
                temp = 3^60
                x = x /  temp
                x = mySin (x)                            ' 通过优化的泰勒公式计算sinx的值  ,步骤见前面。            
                For i = 1 To temp                        '倍角公式返回后运算部分,是不是太简单!
                     x=3x-4x^3
                Next
               







补充内容 (2016-4-18 09:10):
范了一个小错误,FFT 乘法做成DLL时,做成了Debug版,现改成Release,万位乘万位,精度一万用时17ms
Debug版用时37ms,现在万位精度sinx,用时6.5ms

补充内容 (2016-4-19 13:15):
6.5s

补充内容 (2016-4-19 13:27):
for i=1  to   60
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-19 12:13:39 | 显示全部楼层
楼上写错了吧?“万位乘万位,精度一万用时17ms”,“现在万位精度sinx,用时6.5ms”, sin(x)的计算时间怎么能比乘法更快呢?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-19 14:04:00 | 显示全部楼层
你的sin(x)算法仍然不够好,依然有很大的潜力。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-20 15:58:39 | 显示全部楼层
liangbch 发表于 2016-4-19 14:04
你的sin(x)算法仍然不够好,依然有很大的潜力。

你的前面第一个方法我已采用,化简弧度使它小于pi/2,后面两种优化方法已大概看懂,可以叠加在我的方法上,再次加快三角函数运算,只是实现有点复杂,过一段时间看看能不能加上,我的对数函数现在没有头绪,不知道怎么办,你有什么好的建议吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-20 20:04:00 | 显示全部楼层
对于对数函数的任意精度计算,我还是有发言权的。我对对数函数的算法琢磨了将近30年。目前正在探索一种全新的算法并实现一个任意精度对数库,期望其性能超过目前所有的同类数学软件,包括Mathematics,Maple,PARI/GP,MPFP。

我搜索了一下,下面网页留下我的足迹。你可参考下。
http://bbs.pfan.cn/post-241245.html
http://www.cnblogs.com/skyivben/archive/2013/02/15/2912914.html
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-20 20:18:56 | 显示全部楼层
liangbch 发表于 2016-4-20 20:04
对于对数函数的任意精度计算,我还是有发言权的。我对对数函数的算法琢磨了将近30年。目前正在探索一种全新 ...

非常感谢,我一直找不到这方面的资料!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-5-9 18:54:53 | 显示全部楼层
\( \frac {5^4 \times 7} {2 \times 3^7}= \frac{4375}{4374} = 1.000228623685413808870598994 \)

\( \frac{2^{51}}{3^{13} \times 5 \times 7^{10}} =1.000007052909423015206070009 \)
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-10 19:59:32 | 显示全部楼层
本帖最后由 落叶 于 2016-5-10 20:56 编辑

你这个式子竟然不花很大代价获得了很高的自变量分解,采用这个方法提速效果应该很好!正在学习中...
我目前已采用了比较容易编程的方法。就是开方分解(早先调试开方程序时就发现,大于零的数进行开方运算,当开方次数无限增加时,开方值无限接近1),和你的方法可以叠加使用,我的方法如下:
        Ln(x)=Ln(y*10^n)
                =Ln(y)+n*Ln(10)                    'Ln(10)提前算出
        For i = 1 To 45                                '针对一万精度大约进行45次分解
                   y = myKpf(y)                         ‘开平方
        Next
        对y进行泰勒级数运算(没能象sinx一样找到简化算法,现只能老老实实硬算),针对一万精度采用(精度除以13)次泰勒级数运算。
        公式是这个:ln(1+x)=x-x^2/2+x^3/3-……+(-1)^(k-1)*(x^k)/k(|x|
        y = y*(2^45)
        x=y+n*Ln(10)
        myLn = x             '返回值
     经测试万位精度用时32m左右,BTCAL疯狂计算器v2.5用时40m左右,嘿嘿,现在只能以打败它为目标。
     针对你的第二个式子,Ln(9)我需要20次的分解(也就是开1048576次方)用时6s,如果在我的式子上加上你的方法,可以减少6S的运算时间。当要求精度较低时,你的方法提速太恐怖了,可以瞬刷结果,
    不过我还没明白如何具体分解。

补充内容 (2016-5-27 13:53):
万位精度用时32m左右,BTCAL疯狂计算器v2.5用时40m左右
改为:万位精度用时32秒左右,BTCAL疯狂计算器v2.5用时40秒左右
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-29 21:33 , Processed in 0.043044 second(s), 14 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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