C++ で円周率を計算する

「C++で巨大素数を判定・生成する」で使ったGMPライブラリを使用し、チュドノフスキーの公式により、指定した任意桁数の円周率を求めるプログラムを作る。出力量が大きいのでテキスト・ファイルへ行う。メモリがある限りの桁数を計算できるが、実際にメモリ32GバイトのPCを使って1千万桁は余裕で計算できた。
コマンドライン・プログラムなので、Windowsだけではなく、他のOSでも動作するソース・プログラムになっている。

（2024年2月17日）小数点以下15桁未満も計算できるように改良．英語表記追加(-eオプション)，使用ライブラリ等を更新
（2023年10月8日）使用ライブラリをバージョンアップ
（2023年5月28日）使用ライブラリをバージョンアップ

サンプル・プログラム

ダウンロード

圧縮ファイルの内容
pi.exe	実行プログラム本体（CUI版）
pi.cpp	ソース・プログラム
makefile	ビルドファイル

pi.cpp 更新履歴
バージョン	更新日	内容
1.0.2	2023/10/08	使用ライブラリをバージョンアップ
1.0.1	2023/05/28	使用ライブラリをバージョンアップ
1.0.0	2023/03/21	初版

使用ライブラリ

多倍長整数演算するためにライブラリ GMP が必要になる。
コマンドライン・オプションを操作する場合などに、オープンソースのライブラリ Boost C++ライブラリが必要になる。
各々の導入方法等については、「C++ 開発環境の準備」をご覧いただきたい。
コマンドラインオプションの定義、取得については、「解説：コマンドラインオプションの定義と取得 - C++ で巨大素数を判定・生成する」をご覧いただきたい。

ビルド

ビルドするには同梱の "makefile" を利用してほしい。

解説：チュドノフスキーの公式

1914年に、インドの天才数学者ラマヌジャンは円周率を高速に計算する公式を発見した。
この公式を改良し、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

チュドノフスキーの公式をコンピュータ計算で解くには 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 < 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))) \) に落ち着く。

解説：初期値など

   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 APPCMD      "pi"                    //アプリケーション名（CUI版）
  22: #define APPVERSION  "1.0.2"                 //バージョン
  23: #define APPYEAR     "2023"                  //作成年
  24: #define REFERENCE   "https://www.pahoo.org/e-soul/webtech/cpp01/cpp01-02-15.shtm"   // 参考サイト
  25:
  26: #define DEFAULT_NUMBER      100000          //小数点以下の桁数（デフォルト）
  27:
  28: //出力
  29: #define DEFAULT_FILENAME    "pi.txt"        //出力ファイル名（デフォルト）
  30: #define DIGITS_PER_COLUMN   10              //1列の桁数
  31: #define DIGITS_PER_LINE     100             //1行の桁数
  32: #define DIGITS_PER_BLOCK    1000            //1ブロックの桁数
  33:
  34: //エラー・メッセージ格納用
  35: string ErrorMessage;
  36:
  37: //Binary Splitting Methodのための構造体
  38: typedef struct _PQT {
  39:     mpz_class P, Q, T;
  40: } PQT;

初期値は、とくに記載のないかぎり自由に変更できる。
デフォルトの計算桁数や出力ファイル名を定数として用意する。
ファイル出力は、慣例に倣い、DIGITS_PER_COLUMN 桁毎に空白をあけ、DIGITS_PER_LINE 桁毎に改行し、DIGITS_PER_BLOCK 桁毎に空行を追加する。
また、Binary Splitting Methodのための3値構造体 PQT を用意する。

解説：チュドノフスキーの公式の実装

  80: // Chudnovskyクラス =========================================================
  81: class Chudnovsky {
  82:     mpz_class A, B, C, D, E, C3_24;
  83:     unsigned long PREC, N;
  84:     double DIGITS_PER_TERM;
  85: private:
  86:     PQT compPQT(int n1, int n2);
  87: public:
  88:     Chudnovsky();
  89:     double compPi(unsigned long number, string filename);
  90: };
  91:
  92: /*
  93:  * コンストラクタ
  94:  */
  95: Chudnovsky::Chudnovsky() {
  96:     //定数
  97:     A       = 13591409;
  98:     B       = 545140134;
  99:     C       = 640320;
100:     D       = 426880;
101:     E       = 10005;
102:     DIGITS_PER_TERM = 14.1816474627254776555;   //log(53360^3) / log(10)
103:     C3_24   = C * C * C / 24;
104: }
105:
106: /*
107:  * Binary Splitting Method
108:  * @param   int m1, n2
109:  * @return  PQT構造体
110:  */
111: PQT Chudnovsky::compPQT(int n1, int n2) {
112:     int m;
113:     PQT res;
114:
115:     if (n1 + 1 == n2) {
116:         res.P  = (2 * n2 - 1);
117:         res.P *= (6 * n2 - 1);
118:         res.P *= (6 * n2 - 5);
119:         res.Q  = C3_24 * n2 * n2 * n2;
120:         res.T  = (A + B * n2) * res.P;
121:         if ((n2 & 1) == 1) {
122:             res.T = - res.T;
123:         }
124:     } else {
125:         m = (n1 + n2) / 2;
126:         PQT res1 = compPQT(n1, m);
127:         PQT res2 = compPQT(m, n2);
128:         res.P = res1.P * res2.P;
129:         res.Q = res1.Q * res2.Q;
130:         res.T = res1.T * res2.Q + res1.P * res2.T;
131:     }
132:     return res;
133: }
134:
135: /*
136:  * 円周率を計算する．
137:  * @param   unsigned long number 計算する小数点以下桁数
138:  * @param   string filename 円周率を出力するファイル名
139:  * @return  double 計算時間（秒）
140:  */
141: double Chudnovsky::compPi(unsigned long number=DEFAULT_NUMBER, string filename=DEFAULT_FILENAME) {
142:     N       = number / DIGITS_PER_TERM;
143:     PREC    = number * log2(10);
144:
145:     //計算開始時刻
146:     cout << "計算開始... " << number << "桁" << endl;
147:     clock_t t0 = clock();
148:
149:     //円周率を計算する．
150:     PQT PQT = compPQT(0, N);
151:     mpf_class pi(0, PREC);
152:     pi  = D * sqrt((mpf_class)E) * PQT.Q;
153:     pi /= (A * PQT.Q + PQT.T);
154:
155:     //計算完了時刻
156:     clock_t t1 = clock();
157:     cout << "計算完了." << endl;
158:
159:     //ファイル出力
160:     cout << "ファイル出力... " << filename << endl;
161:     ofstream ofs(filename);
162:
163:     //円周率を文字列に変換する．
164:     stringstream ss;
165:     ss.precision(number + 1);
166:     ss << pi;
167:     string str_x = ss.str();
168:
169:     //整数部
170:     ofs << "pi = " << str_x[0] << str_x[1] << endl;
171:     //小数部
172:     for (size_t i = 2; i < str_x.length(); i++) {
173:         ofs << str_x[i];
174:         if ((i - 1) % DIGITS_PER_COLUMN == 0) {
175:             ofs << " ";
176:         }
177:         if ((i - 1) % DIGITS_PER_LINE == 0) {
178:             ofs << endl;
179:         }
180:         if ((i - 1) % DIGITS_PER_BLOCK == 0) {
181:             ofs << endl;
182:         }
183:     }
184:     cout << "ファイル出力完了." << endl;
185:
186:     return (double)(t1 - t0);
187: }
188:
189: // End of class =========================================================

チュドノフスキーの公式と Binary Splitting Method は、前述の定義通りにクラス Chudnovsky に実装した。
計算した円周率はテキスト・ファイルに出力する。デフォルトのファイル名は DEFAULT_FILENAME に設定してあり、オプション -f -file を使って任意のファイル名を指定できる。
ファイルに出力する際、いったん文字列に変換してから、DIGITS_PER_COLUMN 桁毎に空白をあけ、DIGITS_PER_LINE 桁毎に改行し、DIGITS_PER_BLOCK 桁毎に空行を追加する処理を加えている。
なお、チュドノフスキーの公式は小数点以下15桁未満の時には機能しないので、あらかじめ用意した小数点以下15桁の円周率（文字列）から取り出すようにした。

計算時間

このプログラムで1億桁まで計算させた。
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秒かかった。