Last active
November 19, 2018 06:26
-
-
Save komasaru/719d3d789257805979bfad99cdba5bae to your computer and use it in GitHub Desktop.
Fortran 95 source code to compute PI with FMLIB, using Chudnovsky formula.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
!**************************************************** | |
! 円周率の計算 | |
! * 多倍長演算ライブラリ 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