Fortran 95 source code to compute PI with FMLIB, using Chudnovsky formula.

pi_chudnovsky.f95

!****************************************************
! 円周率の計算
! * 多倍長演算ライブラリ FMLIB 使用
! * Chudnovsky の公式使用
!
! date          name            version
! 2018.11.18    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
program pi_chudnovsky
  use fmzm
  implicit none
  ! 各種定数
  integer,      parameter :: SP = kind(1.0)
  integer(SP),  parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  integer(SP),  parameter :: UID   = 10        ! 出力ファイル用 UID
  character(*), parameter :: F_OUT = "PI.txt"  ! 結果出力ファイル
  ! (Digits per term) log(53360 ** 3) / log(10)
  real(DP),     parameter :: DIGIT_T = 14.181647462725476_DP
  ! (for Chudnovsky)
  integer(SP),  parameter :: CHD_A = 13591409
  integer(SP),  parameter :: CHD_B = 545140134
  integer(SP),  parameter :: CHD_C = 640320
  integer(SP),  parameter :: CHD_D = 426880
  integer(SP),  parameter :: CHD_E = 10005
  ! 各種変数
  character(16)    :: fmt_s          ! 結果出力フォーマット
  character(10**9) :: s              ! 結果文字列
  integer(SP)      :: digit          ! 計算桁数（小数点以下）
  integer(SP)      :: ios            ! ファイル IO ステータス
  integer(SP)      :: n              ! N 格納用
  real(DP)         :: t_0, t_1, t_2  ! 実行時間計算用
  type(IM)         :: a, b, c, d, e  ! Chudnovsky 用(A, B, C, D, E)
  type(IM)         :: c3_24          ! C**3 / 24 格納用
  type(FM)         :: pi             ! 円周率
  ! 構造型
  type :: t_pqt
    type(IM) :: p, q, t
  end type t_pqt

  ! 計算桁数取得
  write (*, '(A)', advance="no") "Number of digits to calculate: "
  read (* , '(I12)') digit
  if (digit == 0) stop
  print '(/, A, I0, A)', "#### PI COMPUTATION ( ", digit, "-digit )"

  ! 処理開始
  call FM_SET(digit + 2)

  ! 各種初期値計算
  ! * a, b, c, d, e: Chudnovsky 用
  ! * c3_24: C * C * C / 24 = 10939058860032000
  ! * n: 計算項数
  a = TO_IM(CHD_A)
  b = TO_IM(CHD_B)
  c = TO_IM(CHD_C)
  d = TO_IM(CHD_D)
  e = TO_IM(CHD_E)
  c3_24 = c**3 / 24
  n = int(digit / DIGIT_T)

  ! 実計算
  call cpu_time(t_0)
  call compute(n, pi)
  !call FMPRINT(pi)

  ! 計算終了
  call cpu_time(t_1)
  print '("Duration of computaion: ", F8.3, " sec.")', t_1 - t_0

  ! 書き込み用テキストファイル OPEN
  ! * 以下では、作成済みのファイルを誤消去しないよう status を "new" にしているが、
  !   上書きするようにしたければ、 "replace" にする。
  open (unit   = UID,         &
      & iostat = ios,         &
      & file   = F_OUT,       &
      & action = "write",     &
      & form   = "formatted", &
      & status = "new")
      !& status = "replace")
  if (ios /= 0) then
    print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, F_OUT
    stop
  end if

  ! ファイル出力
  write (fmt_s, '("F", I0, ".", I0)') digit + 2, digit
  call FM_FORM(fmt_s, pi, s)
  write (UID, '(A)') trim(s)

  ! 書き込み用テキストファイル CLOSE
  close(UID)

  ! 処理終了
  call cpu_time(t_2)
  print '("           file output: ", F8.3, " sec.")', t_2 - t_1
  print '("                 total: ", F8.3, " sec.")', t_2 - t_0

  stop
contains
  ! 計算
  ! * 級数展開の計算量削減のために BSA 法(Binary Splitting Algorithm) を使用
  ! * 計算結果を 10 ** digit で除算したものが真の値
  !
  ! :param(in)  integer(4) n: 計算項数
  !: param(out) type(IM)  pi: 円周率
  subroutine compute(n, pi)
    implicit none
    integer(SP), intent(in)  :: n
    type(FM),    intent(out) :: pi
    type(t_pqt) :: pqt

    call bsa(0, n, pqt)
    pi = d * sqrt(TO_FM(e)) * pqt%q / (a * pqt%q + pqt%t)
  end subroutine

  ! BSA 法
  !
  ! :param(in)    integer(4)  n_0: 左端項
  ! :param(in)    integer(4)  n_1: 右端項
  ! :param(inout) type(t_pqt) pqt: 構造型 PQT
  recursive subroutine bsa(n_0, n_1, pqt)
    implicit none
    integer(SP), intent(in)  :: n_0, n_1
    type(t_pqt), intent(out) :: pqt
    type(t_pqt) :: pqt_l, pqt_r  ! 左端〜中間項、中間〜右端項の PQT
    type(IM)    :: in_1          ! n_1 を IM 化したもの
    integer(SP) :: m             ! 中間項

    in_1 = TO_IM(n_1)
    if (n_0 + 1 == n_1) then
      pqt%p = (in_1 * 2 - 1) * (in_1 * 6 - 1) * (in_1 * 6 - 5)
      pqt%q = c3_24 * in_1 * in_1 * in_1
      pqt%t = (a + b * in_1) * pqt%p
      if (iand(n_1, 1) == 1) pqt%t = -pqt%t
    else
      m = int((n_0 + n_1) / 2)
      call bsa(n_0, m, pqt_l)
      call bsa(m, n_1, pqt_r)
      pqt%p = pqt_l%p * pqt_r%p
      pqt%q = pqt_l%q * pqt_r%q
      pqt%t = pqt_l%t * pqt_r%q + pqt_l%p * pqt_r%t
    end if
  end subroutine bsa
end program pi_chudnovsky