几种快速计算任意精度圆周率π的算法

风平老涡

从古至今，计算圆周率π的精度在不断提高，算法也在不断改进。在古代，一般是用园内嵌及园外切多边形来逼近 π。到了近代，通常是用三角函数的泰勒展开来近似，特别是采用反正切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算法

[Fortran] 纯文本查看 复制代码

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算法

[Fortran] 纯文本查看 复制代码

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