Fortran Coder

查看: 63|回复: 0

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

[复制链接]

57

帖子

2

主题

0

精华

专家

F 币
391 元
贡献
170 点

规矩勋章

发表于 2020-9-10 11:14:11 | 显示全部楼层 |阅读模式
本帖最后由 风平老涡 于 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算法
[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



回复

使用道具 举报

您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-9-27 08:50

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表