zeroieme
发表于 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 MOD2pi '弧度化简,x \2pi,取余数,余数为角x,sinx值不变,x 现在小于2pi
x = x MODpi '弧度化简,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=1to 60
liangbch
发表于 2016-4-19 12:13:39
楼上写错了吧?“万位乘万位,精度一万用时17ms”,“现在万位精度sinx,用时6.5ms”, sin(x)的计算时间怎么能比乘法更快呢?
liangbch
发表于 2016-4-19 14:04:00
你的sin(x)算法仍然不够好,依然有很大的潜力。
落叶
发表于 2016-4-20 15:58:39
liangbch 发表于 2016-4-19 14:04
你的sin(x)算法仍然不够好,依然有很大的潜力。
你的前面第一个方法我已采用,化简弧度使它小于pi/2,后面两种优化方法已大概看懂,可以叠加在我的方法上,再次加快三角函数运算,只是实现有点复杂,过一段时间看看能不能加上,我的对数函数现在没有头绪,不知道怎么办,你有什么好的建议吗?
liangbch
发表于 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年。目前正在探索一种全新 ...
非常感谢,我一直找不到这方面的资料!
liangbch
发表于 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秒左右