Fortran source code to compute Pi with arctan's formula.

calc_pi_arctan.f95

!***********************************************************
! 円周率計算 by Arctan 系公式
! ( 各項の Arctan を個別に計算後に加減算する方法 )
!
!    1: Machin
!    2: Klingenstierna
!    3: Eular
!    4: Eular(2)
!    5: Gauss
!    6: Stormer
!    7: Stormer(2)
!    8: Takano
!
! コンパイル方法：
!   $ f95 calc_pi_arctan.f95 -o calc_pi_arctan
!***********************************************************
!
!***********************************************************
! 定数モジュール
!***********************************************************
module constants
  implicit none
  ! 公式名
  character(len=14), parameter :: FORMULA(8) = (/ &
    'Machin        ', 'Klingenstierna', 'Eular         ', 'Eular2        ', &
    'Gauss         ', 'Stormer       ', 'Stormer2      ', 'Takano        ' /)
  ! 公式の項数
  integer, parameter :: KOSU(8) = (/ 2, 3, 2, 3, 3, 3, 4, 4 /)
  ! 公式内係数
  integer, parameter :: KEISU(12, 8) = reshape(source = (/ &
     16, 1,  5,  -4, 1, 239,   0, 0,   0,  0, 0,      0,   &  ! 1: Machin
     32, 1, 10,  -4, 1, 239, -16, 1, 515,  0, 0,      0,   &  ! 2: Klingenstierna
     20, 1,  7,   8, 3,  79,   0, 0,   0,  0, 0,      0,   &  ! 3: Eular
     16, 1,  5,  -4, 1,  70,   4, 1,  99,  0, 0,      0,   &  ! 4: Eular(2)
     48, 1, 18,  32, 1,  57, -20, 1, 239,  0, 0,      0,   &  ! 5: Gauss
     24, 1,  8,   8, 1,  57,   4, 1, 239,  0, 0,      0,   &  ! 6: Stormer
    176, 1, 57,  28, 1, 239, -48, 1, 682, 96, 1,  12943,   &  ! 7: Stormer(2)
     48, 1, 49, 128, 1,  57, -20, 1, 239, 48, 1, 110443    &  ! 8: Takano
  /), shape =(/ 12, 8 /))
  ! 出力ファイル
  character(len=*), parameter :: FILE_PRE = 'PI_'   ! プリフィックス
  character(len=*), parameter :: FILE_EXT = '.txt'  ! 拡張子
  ! 入力メッセージ
  character(len=*), parameter :: STR_INF_1 = "1:Machin, 2:Klingenstierna, 3:Eular, 4:Eular(2),"
  character(len=*), parameter :: STR_INF_2 = "5:Gauss, 6:Stormer, 7:Stormer(2), 8:Takano : "
  character(len=*), parameter :: STR_INF_3 = "Please input number of Pi Decimal-Digits : "
  ! 出力文字列 ( コンソール、テキストファイル共通 )
  character(len=*), parameter :: STR_TITLE   = "** Pi Computation by Arctan method **"
  character(len=*), parameter :: STR_FORMULA = "   Formula = "
  character(len=*), parameter :: STR_DIGITS  = "   Digits  = "
  character(len=*), parameter :: STR_TIME    = "   Time(s) = "
  ! 多桁計算
  integer, parameter :: MAX_DIGITS = 100000000  ! １つの int で８桁扱う
end module constants

!***********************************************************
! 計算モジュール
!***********************************************************
module calc
  use constants
  implicit none
  private              ! モジュール外非公開
  public :: calc_pi    ! calc_pi のみモジュール外公開
  integer :: ary_size  ! 配列サイズ

contains

  !********************
  ! 円周率計算
  !********************
  subroutine calc_pi(f, n)
    integer, intent(in) :: f, n    ! 公式番号、計算桁数
    integer :: nn(KOSU(f))         ! 計算項数
    ! 計算用配列 ( サイズは後に再設定 )
    ! ( 各項の(2k-1)部分を除いた部分、各項の(2k-1)部分含む部分、各項の総和、総和 )
    integer, allocatable :: a(:,:), q(:,:), sk(:,:), s(:)
    integer :: i, k                ! ループインデックス
    integer t1, t2, t_rate, t_max  ! 時間計測
    real :: tt                     ! 総時間

    ! ==== 計算開始時刻
    call system_clock(t1)

    ! ==== 初期設定
    ary_size = (n / 8) + 1  ! 配列サイズ
    do i = 1, KOSU(f)       ! 計算項数
      nn(i) = (n / (log10(real(KEISU((i - 1) * 3 + 3, f))) - &
              log10(real(KEISU((i - 1) * 3 + 2, f)))) + 1) / 2 + 1
    end do

    ! ==== 配列再宣言・初期化
    allocate(a( ary_size + 2, KOSU(f)))  ! 各項の(2k-1)部分を除いた部分
    allocate(q( ary_size + 2, KOSU(f)))  ! 各項の(2k-1)部分含む部分
    allocate(sk(ary_size + 2, KOSU(f)))  ! 各項の総和
    allocate(s( ary_size + 2))           ! 総和
    a = 0; q = 0; sk = 0; s = 0

    ! ==== 計算初期値
    do i = 1, KOSU(f)
      if (KEISU((i - 1) * 3 + 1, f) >= 0) then
        a(1, i) = KEISU((i - 1) * 3 + 1, f) * KEISU((i - 1) * 3 + 3, f)
      else
        a(1, i) = -(KEISU((i - 1) * 3 + 1, f) * KEISU((i - 1) * 3 + 3, f))
      end if
      ! 分子が 1 で無ければ、さらに分子の値で除算しておく
      if (KEISU((i - 1) * 3 + 2, f) /= 1) then
        call long_div(a(:, i), KEISU((i - 1) * 3 + 2, f), a(:, i))
      end if
    end do

    ! ==== 計算
    do i = 1, KOSU(f)
      do k = 1, nn(i)
        ! 各項の分数の部分
        call long_div(a(:, i), KEISU((i - 1) * 3 + 3, f), a(:, i))
        call long_div(a(:, i), KEISU((i - 1) * 3 + 3, f), a(:, i))
        ! 分子が 1 でない時
        if (KEISU((i - 1) * 3 + 2, f) /= 1) then
          call long_mul(a(:, i), KEISU((i - 1) * 3 + 2, f), a(:, i))
          call long_mul(a(:, i), KEISU((i - 1) * 3 + 2, f), a(:, i))
        end if

        ! 各項の 1 / (2 * k - 1) の部分
        call long_div(a(:, i), 2 * k - 1, q(:, i));

        ! 各項の総和に加減算
        if (mod(k, 2) /= 0) then
          call long_add(sk(:, i), q(:, i), sk(:, i))
        else
          call long_sub(sk(:, i), q(:, i), sk(:, i))
        end if
      end do
    end do

    ! ==== 各項の総和の加減算
    do i = 1, KOSU(f)
      if (KEISU((i - 1) * 3 + 1, f) >= 0) then
        call long_add(s, sk(:, i), s)
      else
        call long_sub(s, sk(:, i), s)
      end if
    end do

    ! ==== 計算終了時刻
    call system_clock(t2, t_rate, t_max)

    ! ==== 計算時間
    if ( t2 < t1 ) then
      tt = (t_max - t1 + t2) / dble(t_rate)
    else
      tt = (t2 - t1) / dble(t_rate)
    endif

    ! ==== 結果出力
    call display(f, n, tt, s);
  end subroutine calc_pi

  !********************
  ! ロング + ロング
  !********************
  subroutine long_add(a, b, q)
    integer, intent(in) :: a(:)  ! 引数(被加数配列)
    integer, intent(in) :: b(:)  ! 引数(加数配列)
    integer :: q(:)              ! 引数(計算後配列)(参照渡し)
    integer :: i                 ! ループインデックス
    integer :: cr                ! 繰り上がり

    cr = 0
    do i = ary_size + 2, 1, -1
      q(i) = a(i) + b(i) + cr
      if (a(i) < MAX_DIGITS) then
        cr = 0
      else
         q(i) = a(i) - MAX_DIGITS
         cr = 1
      end if
    end do
  end subroutine long_add

  !********************
  ! ロング - ロング
  !********************
  subroutine long_sub(a, b, q)
    integer, intent(in) :: a(:)  ! 引数(被減数配列)
    integer, intent(in) :: b(:)  ! 引数(減数配列)
    integer :: q(:)              ! 引数(計算後配列)(参照渡し)
    integer :: i                 ! ループインデックス
    integer :: br                ! 借り

    do i = ary_size + 2, 1, -1
      q(i) = a(i) - b(i) - br
      if (a(i) >= 0) then
        br = 0
      else
        q(i) = a(i) + MAX_DIGITS
        br = 1
      end if
    end do
  end subroutine long_sub

  !********************
  ! ロング * ショート
  !********************
  subroutine long_mul(a, b, q)
    integer, intent(in) :: a(:)  ! 引数(被乗数配列)
    integer, intent(in) :: b     ! 引数(乗数)
    integer :: q(:)              ! 引数(計算後配列)(参照渡し)
    integer :: i                 ! ループインデックス
    integer :: cr                ! 繰り上がり
    integer :: w                 ! 被乗数ワーク

    cr = 0
    do i = ary_size + 2, 1, -1
      w = a(i)
      q(i) = mod(w * b + cr, MAX_DIGITS)
      cr = (w * b + cr) / MAX_DIGITS
    end do
  end subroutine long_mul

  !********************
  ! ロング / ショート
  !********************
  subroutine long_div(a, b, q)
    integer, intent(in) :: a(:)  ! 引数(被除数配列)
    integer, intent(in) :: b     ! 引数(除数)
    integer :: q(:)              ! 引数(計算後)(参照渡し)
    integer :: i                 ! ループインデックス
    integer(8) :: r              ! 剰余ワーク
    integer :: w                 ! 被除数ワーク

    r = 0
    do i = 1, ary_size + 2
      w = a(i)
      q(i) = (w + r) / b
      r = mod(w + r, b) * MAX_DIGITS
    end do
  end subroutine long_div

  !********************
  ! 結果出力
  !********************
  subroutine display(f, n, tt, s)
    integer, intent(in) :: f           ! 公式番号
    integer, intent(in) :: n           ! 計算桁数
    real,    intent(in) :: tt          ! 総時間
    integer, intent(in) :: s(:)        ! 計算結果配列
    character(len=14) :: formula_name  ! 公式名
    character(len=21) :: out_file      ! 出力ファイル名
    integer :: i                       ! ループインデックス

    ! ==== 公式名・出力ファイル名設定
    formula_name = FORMULA(f)
    out_file = FILE_PRE // trim(FORMULA(f)) // FILE_EXT

    ! ==== コンソール出力
    print '(/a)',      STR_TITLE
    print '(a,a)',     STR_FORMULA, trim(formula_name)
    print '(a,i10)',   STR_DIGITS,  n
    print '(a,f10.3)', STR_TIME,    tt

    ! ==== ファイル出力
    open (10, file=out_file, status='replace')
    write (10, '(a)')       STR_TITLE
    write (10, '(a,a)')     STR_FORMULA, trim(formula_name)
    write (10, '(a,i10)')   STR_DIGITS,  n
    write (10, '(a,f10.3)') STR_TIME,    tt
    write (10, *)
    write (10, '(i11,".")') s(1)
    do i = 2, ary_size
      if (mod(i, 10) == 2) then
        write (10, '(i8.8,":")', advance='no') (i - 2) * 8 + 1
      end if
      write (10, '(i9.8)', advance='no') s(i)
      if (mod(i, 10) == 1) then
        write (10, *)
      end if
    end do
    close (10)
  end subroutine display
end module calc

!***********************************************************
! 主処理
!***********************************************************
program calc_pi_arctan
  use constants    ! 定数モジュール
  use calc         ! 計算モジュール
  implicit none    ! 暗黙の型指定を使用しない
  integer :: f, n  ! 公式番号、計算桁数

  ! ==== 公式番号入力
  write (*, '(a)')               STR_INF_1
  write (*, '(a)', advance='no') STR_INF_2
  read *, f
  ! ==== 計算桁数入力
  write (*, '(a)', advance='no') STR_INF_3
  read *, n

  ! ==== 円周率計算
  call calc_pi(f, n)
end program calc_pi_arctan