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

落叶

以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方法，（没查到中文资料）
有知道的大大们讲一下，不胜感谢！
有好的算法也贴出来，大家参考参考！

liangbch

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

liangbch

By the way, "万位精度sinx在32位机上用时100秒左右"这个速度太慢了。 我的理想目标是 100万 做到10秒左右。

落叶

想问一下，你的三角函数运行那么快，是因为基础函数速度快，还是对泰勒展开式有什么好的加速方法，
可以简单介绍一下吗？

liangbch

这是我期望的目标，还没有实现，目前只是有些想法而已。因为PiFast 计算百万位 Pi 和百万位e 都可达到10秒以内，所以我期望sin(x)函数也可以达到10秒以内。
我想，即使这个目标实现不了，也不会比这个目标差很多。

liangbch

关于任意精度初等函数快速算法， 理查德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步以后直接运算泰勒公式，具体方法，这里就不再透露了。

liangbch

一点说明：
因为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次大数开平方运算。

落叶

非常感谢你的详细讲解，我在调试三角函数时，也发现了三角函数的取值不同，相同级数运算所获精度不同的问题，网上也有人说三角函数泰勒展开在某些角度逼近太慢，因为一直没查到这
方面的理论知识，一直不知道如何入手，看了你的讲解，我明白了原理（使自变量尽量缩小到0附近，可以极大减少泰勒公式项数，降低了计算量），后面的流程我正在看，里面的步骤，经验
可以极大减少我的盲目摸索，再次感谢！
再请教一下，我的那个方法是以加速泰勒公式为目的，当时想出来时，估计速度提高1-2倍，没想到实现后竟然提高了十倍，不知道能不能先按你所述方法优化后，再按我的方法再次优化加速？

zeroieme

liangbch 发表于 2016-4-13 08:03
化成双阶乘岂不更专业一点，虽然不能加快速度。上式可变为

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

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

liangbch

zeroieme 发表于 2016-4-14 17:24
是不是小除法比一次大除法慢？

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

一楼的算法就是 秦九韶算法。