Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created December 9, 2017 07:28
Embed
What would you like to do?
Fortran source code to compute Pi with arctan's formula.
!***********************************************************
! 円周率計算 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 ! 1つの int で8桁扱う
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment