1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
| !****************************************************
! 円周率の計算
! * 多倍長演算ライブラリ 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) ! <= FMLIB デフォルトのフォーマットで画面出力する場合
! 計算終了
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
|