Chudonovskyの公式

14桁ずつ数値が決まっていきます。

■実行結果

収束の様子と真値との差を表示しています。最後に4*ATAN(1)で真値を書き出しています。

■ソースコード

MODULE m_pi
 IMPLICIT NONE
 INTEGER, PARAMETER :: kq = KIND(1.0q0)

 CONTAINS


  REAL(kq) ELEMENTAL FUNCTION pi_chudonovsky(niter)
   INTEGER, INTENT(IN) :: niter
   REAL(kq) :: a, b, c, d, e, f, s
   INTEGER :: k, m, n

   a = 1.0_kq
   b = 13591409.0_kq
   c = 1.0_kq
   d = 1.0_kq
   e = 1.0_kq
   f = 426880.0_kq * SQRT(10005.0_kq) ! 640320.0_kq**(3.0_kq / 2.0_kq) / 12.0_kq
   m = 0
   n = 0
   s = a * b / (c * d * e)
   DO k = 1, niter
    m = m + 6
    n = n + 3
    a = a * m * (m - 1) * (m - 2) * (m - 3) * (m - 4) * (m - 5)
    b = b + 545140134.0_kq
    c = c * k**3
    d = d * n * (n - 1) * (n - 2)
    e = e * 640320.0_kq**3
    s = s + (-1)**k * a * b / (c * d * e)
   END DO
   pi_chudonovsky = f / s

   RETURN
  END FUNCTION pi_chudonovsky

END MODULE m_pi

!=======================================

PROGRAM test
  USE m_pi
  IMPLICIT NONE
  REAL(kq) :: pi(3)

  pi = pi_chudonovsky([0, 1, 2])
  PRINT '( f40.32)', pi
  PRINT '(es40.32)', pi - 4.0_kq * ATAN(1.0_kq)
  PRINT *
  PRINT '(f40.32)', 4.0_kq * ATAN(1.0_kq)
  
  STOP
END PROGRAM test