Fortran Coder

标题: 几种快速计算任意精度圆周率π的算法 [打印本页]

作者: 风平老涡 时间: 2020-9-10 11:14
标题: 几种快速计算任意精度圆周率π的算法
本帖最后由风平老涡于 2020-9-10 21:16 编辑

从古至今，计算圆周率π的精度在不断提高，算法也在不断改进。在古代，一般是用园内嵌及园外切多边形来逼近 π。到了近代，通常是用三角函数的泰勒展开来近似，特别是采用反正切arctan函数。这类采用arctan泰勒展开的方法又称为Machin法，其中最简单的公式如下：
[latex]\frac{\pi }{4} = 4 arctan\frac{1}{5} - arctan\frac{1}{239}[/latex]
请参考源代码（http://piany.w.fcode.cn/）。

而到了现代（1988年），则产生了号称最快的Chudnovsky兄弟算法。这个算法是基于印度神人Ramanujan的π公式，其表达式如下：
[latex]\frac{1}{\pi } = 12 \sum_{q=0}^{\infty }\frac{(-1)^{q}(6q)!(545140134q+13591409)}{(3q)!(q!)^{3}(640320)^{3q+3/2}}[/latex]

以上这些算法都是线性收敛。为了达到快速计算高精度的圆周率，又产生了高阶收敛的迭代公式，如二阶收敛的Gauss-Legendre算法, 四阶收敛的Borwein算法等。以下为Borwein计算公式：
[latex]a_{0} = 2 (\sqrt{2} - 1)^{2}[/latex]
[latex]y_{0} = \sqrt{2} - 1[/latex]
[latex]y_{k+1} = \frac{1 - (1 - y_{k}^{4})^{1/4}}{1 + (1 - y_{k}^{4})^{1/4}}[/latex]
[latex]a_{k+1} = a_{k} (1 + y_{k+1})^{4} - 2^{2k+3} y_{k+1} (1+y_{k+1} + y_{k+1}^2)[/latex]

自从有了任意精度数据库（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

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)