目次
サンプル・プログラム
| pi.exe | 実行プログラム本体（CUI版） |
| pi.cpp | ソース・プログラム |
| makefile | ビルドファイル |
| バージョン | 更新日 | 内容 |
|---|---|---|
| 1.2.3 | 2025/02/22 | 使用ライブラリ更新 |
| 1.2.2 | 2024/11/03 | 使用ライブラリ更新 |
| 1.2.1 | 2024/07/07 | 使用ライブラリ更新 |
| 1.2.0 | 2024/02/17 | 小数点以下15桁未満も計算できるように改良 |
| 1.1.0 | 2024/02/17 | 英語表記追加(-eオプション)，使用ライブラリ等を更新 |
使用ライブラリ
コマンドライン・オプションを操作する場合などに、オープンソースのライブラリ Boost C++ライブラリが必要になる。
各々の導入方法等については、「C++ 開発環境の準備」をご覧いただきたい。
コマンドラインオプションの定義、取得については、「解説：コマンドラインオプションの定義と取得 - C++ で巨大素数を判定・生成する」をご覧いただきたい。
ビルド
解説：チュドノフスキーの公式
この公式を改良し、1988年にチュドノフスキー兄弟が発表したのが「チュドノフスキー（Chudnovsky）の公式」である。この公式を用いて、1996年に当時の世界記録である80億桁を達成し、2022年3月21日には100兆桁の円周率が計算されている。
\[
\displaystyle
\begin{gather}
\frac{1}{\pi}=\frac{12}{\sqrt{c^3}} \sum_{k=0}^\infty \frac{(-1)^k(6k)!(A + Bk)}{(3k)!(k!)^3C^{3k}} = \frac{12}{\sqrt{c^3}} \sum_{k=0}^\infty P_k
\\
\begin{cases}
A = 13591409 \\
B = 545140134 \\
C = 640320
\end{cases}
\end{gather}
\]
現実には∞まで計算できないのだが、1項目進むごとに
\[ \displaystyle \frac{P_{k+1}}{P_k} = \frac{-24(6k+1)(2k+1)(6k+5)(A+Bk+B)}{C^3(k+1)^3(A+Bk)} \]
だけ小さくなっていく。ここで k→∞ として計算すると
\[ \displaystyle \lim_{k \to \infty} \frac{P_{k+1}}{P_k} = \frac{-24 \times 6^2 \times 2}{C^3} = \frac{-1}{151931373056000} \]
ここで
\[ \displaystyle log_{10} 151931373056000 = 14.1816474627... \]
なので、1項目進むごとに約14桁ずつ精度が高くなっていくことになる。
解説：Binary Splitting Method
チュドノフスキーの公式は無限級数のため、初項から順に足し合わせていくと莫大な時間を要する。そこで、おおむね2分割してトーナメント方式で足し合わせていこうというのが Binary Splitting Method である。
概念としては、
\[
\displaystyle
\frac{P(0,N)}{Q(0,N)} = \sum_{n=0}^{N-1} \frac{a(n)p(0)...p(n)}{b(n)q(0)...q(n)}
\]
の形になるよう、適切な式 \( P(l,u) \), \( Q(l,u) \), \( R(l,u) \) を設定すると、\( l \le m \lt u \) となる \( m \) を用いて、次のように再帰計算ができる。 \( l+1=u \) となる場合は、定義式に \( n=u \) を代入する。
\[
\displaystyle
\begin{eqnarray}
P(l,u) &=& P(l,m)Q(m,u)+R(l,m)P(m,u) \\
Q(l,u) &=& Q(l,m)Q(m,u) \\
R(l,u) &=& R(l,m)R(m,u)
\end{eqnarray}
\]
このとき \( \displaystyle m \simeq \lfloor \frac{l+u}{2} \rfloor \) とすることで、計算量は \( O(log(n) \times M(n\ log(n))) \) に落ち着く。
解説：初期値など
pi.cpp
9: // 初期値 ============================================================
10: #include <iostream>
11: #include <fstream>
12: #include <gmpxx.h>
13: #include <boost/program_options.hpp>
14: #include <boost/format.hpp>
15:
16: using namespace std;
17: using namespace boost::program_options;
18:
19: #define MAKER "pahoo.org" // 作成者
20: #define APPNAMEJP "円周率を求める" // アプリケーション名
21: #define APPNAMEEN "Calculating Pi" // アプリケーション名（英語）
22: #define APPCMD "pi" // アプリケーション名（CUI版）
23: #define APPVERSION "1.2.3" // バージョン
24: #define APPYEAR "2023-2025" // 作成年
25: #define REFERENCE "https://www.pahoo.org/e-soul/webtech/cpp01/cpp01-02-15.shtm" // 参考サイト
26:
27: #define DEFAULT_NUMBER 100000 // 小数点以下の桁数（デフォルト）
28:
29: // 出力ファイルの書式
30: #define DEFAULT_FILENAME "pi.txt" // 出力ファイル名（デフォルト）
31: #define DIGITS_PER_COLUMN 10 // 1列の桁数
32: #define DIGITS_PER_LINE 100 // 1行の桁数
33: #define DIGITS_PER_BLOCK 1000 // 1ブロックの桁数
34:
35: // エラー・メッセージ格納用
36: string ErrorMessage;
37:
38: // Binary Splitting Methodのための構造体
39: typedef struct _PQT {
40: mpz_class P, Q, T;
41: } PQT;
デフォルトの計算桁数や出力ファイル名を定数として用意する。
ファイル出力は、慣例に倣い、DIGITS_PER_COLUMN 桁毎に空白をあけ、DIGITS_PER_LINE 桁毎に改行し、DIGITS_PER_BLOCK 桁毎に空行を追加する。
また、Binary Splitting Methodのための3値構造体 PQT を用意する。
解説：チュドノフスキーの公式の実装
pi.cpp
119: // Chudnovskyクラス =========================================================
120: class Chudnovsky {
121: mpz_class A, B, C, D, E, C3_24;
122: unsigned long PREC, N;
123: double DIGITS_PER_TERM;
124: private:
125: PQT compPQT(int n1, int n2);
126: public:
127: Chudnovsky();
128: double compPi(unsigned long number, string filename, bool en);
129: };
130:
131: /*
132: * コンストラクタ
133: */
134: Chudnovsky::Chudnovsky() {
135: // 定数
136: A = 13591409;
137: B = 545140134;
138: C = 640320;
139: D = 426880;
140: E = 10005;
141: DIGITS_PER_TERM = 14.1816474627254776555; // log(53360^3) / log(10)
142: C3_24 = C * C * C / 24;
143: }
144:
145: /*
146: * Binary Splitting Method
147: * @param int m1, n2
148: * @return PQT構造体
149: */
150: PQT Chudnovsky::compPQT(int n1, int n2) {
151: int m;
152: PQT res;
153:
154: if (n1 + 1 == n2) {
155: res.P = (2 * n2 - 1);
156: res.P *= (6 * n2 - 1);
157: res.P *= (6 * n2 - 5);
158: res.Q = C3_24 * n2 * n2 * n2;
159: res.T = (A + B * n2) * res.P;
160: if ((n2 & 1) == 1) {
161: res.T = - res.T;
162: }
163: } else {
164: m = (n1 + n2) / 2;
165: PQT res1 = compPQT(n1, m);
166: PQT res2 = compPQT(m, n2);
167: res.P = res1.P * res2.P;
168: res.Q = res1.Q * res2.Q;
169: res.T = res1.T * res2.Q + res1.P * res2.T;
170: }
171: return res;
172: }
173:
174: /*
175: * 円周率を計算する．
176: * @param unsigned long number 計算する小数点以下桁数
177: * @param string filename 円周率を出力するファイル名
178: * @param bool en=TRUE：英語表示，false：日本語表示（省略時）
179: * @return double 計算時間（秒）
180: */
181: double Chudnovsky::compPi(unsigned long number=DEFAULT_NUMBER, string filename=DEFAULT_FILENAME, bool en=false) {
182: N = number / DIGITS_PER_TERM;
183: PREC = number * log2(10);
184: string str_x;
185:
186: // 計算開始時刻
187: if (! en) {
188: cout << "計算開始... " << number << "桁" << endl;
189: } else {
190: cout << "Calculation start... " << number << "digits" << endl;
191: }
192: clock_t t0 = clock();
193:
194: // 円周率を計算する．
195: if (number >= 15) {
196: PQT PQT = compPQT(0, N);
197: mpf_class pi(0, PREC);
198: pi = D * sqrt((mpf_class)E) * PQT.Q;
199: pi /= (A * PQT.Q + PQT.T);
200: // 円周率を文字列に変換する．
201: stringstream ss;
202: ss.precision(number + 1);
203: ss << pi;
204: str_x = ss.str();
205: // 小数点以下15桁未満の場合
206: } else {
207: const string piStr = "141592653589793";
208: str_x = "3." + piStr.substr(0, number);
209: }
210:
211: // 計算完了時刻
212: clock_t t1 = clock();
213: if (! en) {
214: cout << "計算完了." << endl;
215: } else {
216: cout << "Calculation complete" << endl;
217: }
218:
219: // ファイル出力
220: if (! en) {
221: cout << "ファイル出力... " << filename << endl;
222: } else {
223: cout << "File output... " << filename << endl;
224: }
225: ofstream ofs(filename);
226:
227:
228: // 整数部
229: ofs << "pi = " << str_x[0] << str_x[1] << endl;
230: // 小数部
231: for (size_t i = 2; i < str_x.length(); i++) {
232: ofs << str_x[i];
233: if ((i - 1) % DIGITS_PER_COLUMN == 0) {
234: ofs << " ";
235: }
236: if ((i - 1) % DIGITS_PER_LINE == 0) {
237: ofs << endl;
238: }
239: if ((i - 1) % DIGITS_PER_BLOCK == 0) {
240: ofs << endl;
241: }
242: }
243: if (! en) {
244: cout << "ファイル出力完了." << endl;
245: } else {
246: cout << "File output completed." << endl;
247: }
248:
249: return (double)(t1 - t0);
250: }
251:
252: // End of class =========================================================
計算した円周率はテキスト・ファイルに出力する。デフォルトのファイル名は DEFAULT_FILENAME に設定してあり、オプション -f -file を使って任意のファイル名を指定できる。
ファイルに出力する際、いったん文字列に変換してから、DIGITS_PER_COLUMN 桁毎に空白をあけ、DIGITS_PER_LINE 桁毎に改行し、DIGITS_PER_BLOCK 桁毎に空行を追加する処理を加えている。
なお、チュドノフスキーの公式は小数点以下15桁未満の時には機能しないので、あらかじめ用意した小数点以下15桁の円周率（文字列）から取り出すようにした。
計算時間
9年前の iMac（Core i5-4670，4コア, 3.4GHz）では132秒、5年前の Let's note CF-SZ6（Core i5-7300U，2コア, 2.6GHz）では213秒、最新の VersaPro UltraLite タイプVN（Core i5-1235U，10コア）では116秒かかった。
参考サイト
- 円周率.jp
- C++ 開発環境の準備：ぱふぅ家のホームページ
コマンドライン・プログラムなので、Windowsだけではなく、他のOSでも動作するソース・プログラムになっている。
（2025年2月22日）使用ライブラリ更新
（2024年11月3日）使用ライブラリ更新
（2024年7月7日）使用ライブラリ更新