数学研发论坛

 找回密码
 欢迎注册
查看: 2687|回复: 36

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

[复制链接]
发表于 2016-4-12 22:22:05 | 显示全部楼层 |阅读模式

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

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

x
以sin()函数为例谈谈我的两种实现:公式是
`\D\sin x= x- \frac{x^3}{3!}+ \frac{x^5}{5!} -\frac{ x^7}{7!}+ \frac{x^9}{9!} . =\sum_{n=1}^5(-1)^n \frac{x^{(2n-1)}}{(2n-1)!}`
第一种是在硬算的基础上简单优化:下面是流程码
               x2 = x*x                                                                     'x*x
               temp1 = x                                                                   'x
               n1=2
               n2=3
               jc=1
               sum=x
               fhbz = 0                                                                       '作为判断加减操作的标志
               Do
                    temp1 = x2*temp1                                                 'x^3,      x^5...
                    jc =jc*n1*n2                                                          '1*2*3,      n!
                    temp3 = temp1/ jc                                                  'x^3/3!,    x^(2n+1)/(2n+1)!
                    If fhbz = 0 Then
                        sum = sum- temp3                                             'x-x^3/3!
                       fhbz = 1
                    Else
                        sum = sum+temp3                                              'x + x^5/5!
                        fhbz = 0
                    End If
                    n1 = n1 + 2
                    n2 =  n2  + 2
                Loop Until n2 = jd Or n1 = jd
因为速度太慢,又改进了一下
化成`\sin x= x- x\cdot x^2(4·5·6·7·8·9 - x^2(6·7·8·9 - x^2(8·9- x^2)))/9!` ,速度提高十倍
流程如下:
                X2 = x*x
                xjc = n!
                n12 = 1
                n1=n
                n2=n-1
                xx = X2
                For i = jd To 5 Step -2                                                 'jd是精度控制参数,这里是n
                    n12 = n1* n2*n12
                    n1 = n1 - 2
                    n2 = n2 - 2
                    xx = n12-xx
                    xx = X2* xx
                Next
                xx = x * xx
                xx = xx / xjc
                sum = x - xx
第二种方法用的高精加,减,乘(万进制硬乘配千位FFT乘),除(估商配迭代除法)是我自已写的,万位精度sinx在32位机上用时100秒左右,
网上有个(BTCAL疯狂计算器v2.5,里面很多运算非常强大)万位精度sinx在我的电脑上用时22秒,不过它的高精乘法速度是我的乘法的十倍,也就是说它的三角函数算法速度可能和我的差不多。

有人建议我用binary spliting方法加速泰勒展开式运算,但是我并不懂什么是binary spliting方法,(没查到中文资料)
有知道的大大们讲一下,不胜感谢!
有好的算法也贴出来,大家参考参考!







毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-13 08:03:18 | 显示全部楼层
sinx= x- x*x^2(4*5*6*7*8*9 - x^2(6*7*8*9 - x^2(8*9- x^2)))/9!

化成双阶乘岂不更专业一点,虽然不能加快速度。上式可变为
sinx= x- x*x^2(5*7*9 - x^2(7*9 - x^2(9- x^2)))/9!!


注: 若n为奇数 n!!= 1*3*5*7 * ...n
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-13 08:04:46 | 显示全部楼层
By the way, "万位精度sinx在32位机上用时100秒左右"这个速度太慢了。 我的理想目标是 100万 做到10秒左右。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-13 08:35:52 | 显示全部楼层
想问一下,你的三角函数运行那么快,是因为基础函数速度快,还是对泰勒展开式有什么好的加速方法,
可以简单介绍一下吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-13 14:12:09 | 显示全部楼层
这是我期望的目标,还没有实现,目前只是有些想法而已。因为PiFast 计算百万位 Pi 和百万位e 都可达到10秒以内,所以我期望sin(x)函数也可以达到10秒以内。
我想,即使这个目标实现不了,也不会比这个目标差很多。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-13 16:21:09 | 显示全部楼层
关于任意精度初等函数快速算法, 理查德P.布伦特(Richard P. Brent)是这方面的专家。关于具体细节,有可查阅下面2篇论文。在论文1中,作者说
若M(p)表示2个p比特数的乘法运算(结果精确到p比特)的复杂度,那么,当n趋于无穷大时,得到n位结果的初等函数(exp,log,argtan,sin,conh等)的算法的复杂度O(M(n)log(n)).
1.Fast multiple-precision evaluation of elementary functions
2.Multiple-precision zero-finding methods and the complexity of elementary function evaluation

下面,我们不谈这位科学家的算法,说说我自己的算法。
1.对于任意实数x,我们很容易求得x0,使得0<x0<Pi/2, 那么sin(x)=sin(x0) 或者sin(x)=-sin(x0),这是初中阶段的内容,不解释。
2.对于0<x0<Pi/2,我们可以找到一个x1, 使得x0=(Pi/60*k)+x1, x1<Pi/60,k是15以内的正整数. Pi/60等于3度,所有3度倍数的正余弦函数值可以精确求出,参见本站帖子 http://bbs.emath.ac.cn/thread-4021-1-1.html
3.现在x1的值总是小于0.052359877

4.计算出以下4个值,这4个值很容易算出而且计算速度非常快。那么,sin(y1)=0.002,sin(y2)=0.00002。这几个值可以存储在文件中,在程序启动的时候载入,也可以在程序启动的实时计算。
  y1=arcsin(0.002)
  y2=arcsin(0.00002)
  y3=arcsin(0.0000002)
  y4=arcsin(0.000000002)

5.将x1拆分为 x1=(k1*y1)+(k2*y2)+(k3*y3)+(k4*y4)+x2, k1,k2,k3,k4为99以内的奇数。
  当k为奇数是,sin(ky)可以表示为sin(y)的多项式,且多项式的系数为整数。参见http://mathworld.wolfram.com/Multiple-AngleFormulas.html
  我们可将sin(ky)的sin(y)多项式的系数预先计算出来,并存储在程序中。 从3到99倍角共49个公式。
  
  按如下顺序计算,共需12次大数乘法,4次大数平方,4次大数开平方
  1. u=sin(k1*y1),//使用多倍角公式,工作量可忽略,下同
  2. v=sqrt(1-u*u)
  3. s1=u, c1=v
  
  4. u=sin(k2*y2),
  5. v=sqrt(1-u*u)
  6. s2=s1*v+c1*u,  //sin(α+β)=sinαcosβ+cosαsinβ
  7. c2=c1*v-s1*u,  //cos(α+β)=cosαcosβ-sinαsinβ
  
  8. u=sin(k3*y3),
  9. v=sqrt(1-u*u)
  10.s3=s2*v+c2*u,  
  11.c3=c2*v-s2*u,  
  
  12.u=sin(k4*y4),
  13.v=sqrt(1-u*u)
  14.s4=s3*v+c3*u
  15.c4=c3*v-s3*u

6.令z=(k1*y1)+(k2*y2)+(k3*y3)+(k4*y4),
   此时sin(x1)=sin(z+x2)=sin(z)cos(x2)+cos(z)sin(x2)=s4*cos(x2)+c4*sin(x2), 而x2小于y4;

通过1-6步,我们将问题sin(x)转化为计算sin(x2),原始的x是任意实数,而x2<y4,此时直接应用泰勒公式,其收敛的速度将非常快。

7.By the way,我的算法将继续缩小自变量的范围,而不是在第6步以后直接运算泰勒公式,具体方法,这里就不再透露了。

评分

参与人数 1威望 +2 金币 +2 贡献 +2 鲜花 +2 收起 理由
落叶 + 2 + 2 + 2 + 2 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-13 17:04:53 | 显示全部楼层
一点说明:
因为u是有限精度数,其精度并不会随着目标精度的而增加。比如当计算精度从1万位扩展到10万时,u的精度不会随之增加,因此
1.在计算“v=sqrt(1-u*u)”,其平方运算的复杂度低于普通的大数平方运算。
2.同理,在计算"s3=s2*v+c2*u",和."c3=c2*v-s2*u",其复杂度也低于普通的大数乘法。
3.因为u是很小的数,故sqrt(1-u*u)有快速算法。
故,最终复杂度实际上是介于6-12次大数乘法,0-4次大数平方,0-4次大数开平方运算。

评分

参与人数 1威望 +2 金币 +2 贡献 +2 鲜花 +2 收起 理由
落叶 + 2 + 2 + 2 + 2 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-4-13 20:40:18 | 显示全部楼层
本帖最后由 落叶 于 2016-4-13 23:19 编辑

非常感谢你的详细讲解,我在调试三角函数时,也发现了三角函数的取值不同,相同级数运算所获精度不同的问题,网上也有人说三角函数泰勒展开在某些角度逼近太慢,因为一直没查到这
方面的理论知识,一直不知道如何入手,看了你的讲解,我明白了原理(使自变量尽量缩小到0附近,可以极大减少泰勒公式项数,降低了计算量),后面的流程我正在看,里面的步骤,经验
可以极大减少我的盲目摸索,再次感谢!
再请教一下,我的那个方法是以加速泰勒公式为目的,当时想出来时,估计速度提高1-2倍,没想到实现后竟然提高了十倍,不知道能不能先按你所述方法优化后,再按我的方法再次优化加速?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-14 17:24:45 | 显示全部楼层
本帖最后由 zeroieme 于 2016-4-14 17:39 编辑
liangbch 发表于 2016-4-13 08:03
化成双阶乘岂不更专业一点,虽然不能加快速度。上式可变为


是不是小除法比一次大除法慢?

x*(1-(x^2/6)*(1-(x^2/20)*(1-(x^2/42)*(1-x^2/72*(……)))))
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-4-14 18:23:00 | 显示全部楼层
zeroieme 发表于 2016-4-14 17:24
是不是小除法比一次大除法慢?

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

一楼的算法就是 秦九韶算法。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2018-1-18 09:58 , Processed in 0.079584 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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