|
|
发表于 2023-10-2 16:56:18
|
显示全部楼层
关于任意精度初等函数快速算法， 理查德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步以后直接运算泰勒公式，具体方法，这里就不再透露了。 |
|
