Borwein Pi ４次公式

参照：http://crd.lbl.gov/~dhbailey/dhbpapers/pi-quest.pdf

■実行結果

■ソースコード

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

 CONTAINS

 ELEMENTAL REAL(kq) FUNCTION pi_Borwein4(n) !P.Borwein : quartic
   INTEGER, INTENT(IN) :: n
   REAL(kq) :: a, y, t, f
   INTEGER :: i

   y = SQRT(2.0_kq) - 1.0_kq
   a = 6.0_kq - 4.0_kq * SQRT(2.0_kq)
   f = 2.0_kq
   DO i = 1, n
     t = SQRT( SQRT(1.0_kq - y**4) )
     y = (1.0_kq - t) / (1.0_kq + t)
     t = (1.0_kq + y)**2
     f =  4.0_kq * f
     a = a * t**2 - f * y * (t - y) 
   END DO
   pi_Borwein4 = 1.0_kq / a

   RETURN
 END FUNCTION pi_Borwein4
!

END MODULE m_pi

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

PROGRAM pi
 USE m_pi
 IMPLICIT NONE

 PRINT *, pi_Borwein4([1, 2, 3])

 STOP
END PROGRAM pi