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

zeroieme

liangbch 发表于 2016-4-14 18:23
一楼的算法就是 秦九韶算法。

差别是除以阶乘，是每次小乘法最后一次大除法还是每次小除法。

落叶

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的增量,我想把它分解出来,看看能不能简化编程,结果最后和秦九韶算法效果一样,这只是一个巧合,不过秦九韶算法真的不错!

落叶

经过几天的学习，查资料，对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

liangbch

楼上写错了吧？“万位乘万位，精度一万用时17ms”，“现在万位精度sinx，用时6.5ms”， sin(x)的计算时间怎么能比乘法更快呢？

liangbch

你的sin(x)算法仍然不够好，依然有很大的潜力。

落叶

liangbch 发表于 2016-4-19 14:04
你的sin(x)算法仍然不够好，依然有很大的潜力。

你的前面第一个方法我已采用,化简弧度使它小于pi/2,后面两种优化方法已大概看懂,可以叠加在我的方法上,再次加快三角函数运算,只是实现有点复杂,过一段时间看看能不能加上,我的对数函数现在没有头绪，不知道怎么办，你有什么好的建议吗？

liangbch

对于对数函数的任意精度计算，我还是有发言权的。我对对数函数的算法琢磨了将近30年。目前正在探索一种全新的算法并实现一个任意精度对数库，期望其性能超过目前所有的同类数学软件，包括Mathematics，Maple，PARI/GP,MPFP。

我搜索了一下，下面网页留下我的足迹。你可参考下。
http://bbs.pfan.cn/post-241245.html
http://www.cnblogs.com/skyivben/archive/2013/02/15/2912914.html

落叶

liangbch 发表于 2016-4-20 20:04
对于对数函数的任意精度计算，我还是有发言权的。我对对数函数的算法琢磨了将近30年。目前正在探索一种全新 ...

非常感谢，我一直找不到这方面的资料！

liangbch

\( \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 \)

落叶

你这个式子竟然不花很大代价获得了很高的自变量分解，采用这个方法提速效果应该很好！正在学习中...
我目前已采用了比较容易编程的方法。就是开方分解(早先调试开方程序时就发现，大于零的数进行开方运算，当开方次数无限增加时，开方值无限接近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秒左右