Skip to content

Instantly share code, notes, and snippets.

@komasaru komasaru/pi_chudnovsky.f95
Last active Nov 19, 2018

Embed
What would you like to do?
Fortran 95 source code to compute PI with FMLIB, using Chudnovsky formula.
!****************************************************
! 円周率の計算
! * 多倍長演算ライブラリ 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.