几种快速计算任意精度圆周率π的算法
本帖最后由 风平老涡 于 2020-9-10 21:16 编辑从古至今，计算圆周率π的精度在不断提高，算法也在不断改进。在古代，一般是用园内嵌及园外切多边形来逼近 π。到了近代，通常是用三角函数的泰勒展开来近似，特别是采用反正切arctan函数。这类采用arctan泰勒展开的方法又称为Machin法，其中最简单的公式如下：
\frac{\pi }{4} = 4 arctan\frac{1}{5} - arctan\frac{1}{239}
请参考源代码（http://piany.w.fcode.cn/）。
而到了现代（1988年），则产生了号称最快的Chudnovsky兄弟算法。这个算法是基于印度神人Ramanujan的π公式，其表达式如下：
\frac{1}{\pi } = 12 \sum_{q=0}^{\infty }\frac{(-1)^{q}(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}(640320)^{3q+3/2}}
以上这些算法都是线性收敛。为了达到快速计算高精度的圆周率，又产生了高阶收敛的迭代公式，如二阶收敛的Gauss-Legendre算法, 四阶收敛的Borwein算法等。以下为Borwein计算公式：
a_{0} = 2 (\sqrt{2} - 1)^{2}
y_{0} = \sqrt{2} - 1
y_{k+1} = \frac{1 - (1 - y_{k}^{4})^{1/4}}{1 + (1 - y_{k}^{4})^{1/4}}
a_{k+1} = a_{k} (1 + y_{k+1})^{4} - 2^{2k+3} y_{k+1} (1+y_{k+1} + y_{k+1}^2)
自从有了任意精度数据库（http://bbs.fcode.cn/thread-2494-1-1.html），计算任意精度圆周率π就变得很容易了。这里给出了Chudnovsky算法及Borwein算法的源代码和计算结果。
1. Chudnovsky算法
program pi_chudnovsky
use fmzm
implicit none
type(im) :: L, L0, X, X0, K
type(fm) :: M, C, t, P1, P2, error
integer(kind=4) :: i, n
character(:), allocatable :: str
character(20) :: fmt, tmp1, tmp2
print *, "Input required number of digits = "
read(*,*) n
write(tmp1, *) n + 3
write(tmp2, *) n - 3
write(fmt, *) "F"//trim(adjustl(tmp1))//"."//trim(adjustl(tmp2))
call fm_set(n)
allocate(character(n) :: str)
L = to_im('13591409')
L0 = to_im('545140134')
X = 1
X0 = to_im('-262537412640768000')
K = 6
M = to_fm('1')
C = to_fm('426880') * sqrt(to_fm('10005'))
t = to_fm(L)
P1 = C / t
error = to_fm(10)**(-n)
i = 0
do
i = i + 1
M = M * to_fm((K * K - 16 ) * K) / to_fm(i)**3
K = K + 12
L = L + L0
X = X * X0
t = t + M * L / X
P2 = C / t
if(abs(P2 - P1) < error) exit
P1 = P2
end do
call fm_form(fmt, P2, str)
write(*,*) "Pi= ", str
end program pi_chudnovsky
2. Borwein算法
program pi_borwein
use fmzm
implicit none
type(fm) :: y, a, t, P1, P2, one, two, four, error
integer :: i, n
character(:), allocatable :: str
character(50) :: fmt, tmp1, tmp2
print *, "Input required number of digits = "
read(*,*) n
write(tmp1, *) n + 3
write(tmp2, *) n - 3
write(fmt, *) "F"//trim(adjustl(tmp1))//"."//trim(adjustl(tmp2))
call fm_set(n)
allocate(character(n) :: str)
one = to_fm(1)
two = to_fm(2)
four = to_fm(4)
y = sqrt(two) - one
a = 6 - four * sqrt(two)
P1 = one / a
error = to_fm(10)**(-n)
i = 0
do
i = i + 1
t = sqrt(one - y ** 4)
t = sqrt(t)
y = (one - t) / (one + t)
t = one + y
a = a * t**4 - two * y * four ** i * (one + y * t)
P2 = one / a
if(abs(P2 - P1) < error) exit
P1 = P2
end do
call fm_form(fmt, P2, str)
write(*,*) "Pi= ", str
end program pi_borwein
结果如下：
(Intel(R) Core(TM) i7-4710HQ CPU @ 2.50GHz, GNU GFortran 10.1.0)
π精度 Chudnovsky Borwein
（位数字） （时间秒） （时间秒）
3000 0.13 0.02
30000 15.23 0.18
300000 2230 2.79
3000000 - 76.98
页:
[1]