Rust による円周率計算プログラム
著者：関 勝寿
公開日：2024年3月17日 - 最終更新日：2024年7月2日
キーワード：
math
rust
python
円周率πを任意桁数計算するプログラム compute-pi をRustで作成した。16GBのメモリを持つPCでは、円周率を3億2000万桁まで計算できる。
計算方法
ガウス・ルジャンドルのアルゴリズム、すなわち初期値を
\[a_0 =1, b_0 =\frac{1}{\sqrt{2}}, t_0 =\frac{1}{4}, p_0 =1\]反復式を
\[\begin{align} a_{n+1} &=\frac{a_n +b_n}{2} \\ b_{n+1} &=\sqrt{a_n b_n} \\ t_{n+1} &=t_n -p_n(a_n - a_{n+1})^2 \\ p_{n+1} &= 2p_n \end{align}\]とすると、円周率πは
\[\pi \approx \frac{(a+b)^2}{4t}\]と近似される、というアルゴリズムにより円周率を任意桁数計算する。多倍長計算のrug crateを使った。
MacBook (Apple M1, 16 GB) では、100万桁を1.5秒で、3億2000万桁を24分で計算できて、これがこのプログラムとマシン能力の限界桁数となる。このプログラムによる円周率3億2000万桁の計算結果は、y-cruncher と Chudnovsky アルゴリズムによる計算結果と完全に一致することが確認された。
高橋大介と金田康正は、1998年にガウス・ルジャンドル法によって分散メモリ型並列計算機による円周率の515億桁計算をした。さらに高橋は2009年にガウス・ルジャンドル法で2兆5769億8037万桁を計算した。以下に、いくつかの参考サイトを示す。
- 世界記録は31兆桁！ 日本人も活躍する円周率「π」計算の最先端 - 柳谷晃, 2021/7/17
- 円周率を求める公式・アルゴリズム - 円周率.jp
- Chudnovsky の公式を用いた円周率の計算用メモ - Peria Peria, 2015/12/9
- パソコンによる円周率小数点以下５兆桁の計算 - 近藤茂, 2011/6/30
- Even more pi in the sky: Calculating 100 trillion digits of pi on Google Cloud - Emma Haruka Iwao, 2022/6/9
- 円周率の歴史 - Wikipedia
使い方
Rustをインストールしてcargoが使えるようになったら、次のようにcompute-pi crateをインストールする。
cargo install compute-pi
すると
compute-pi <桁数>
で、その桁数の円周率が表示される。たとえば、
compute-pi 10
とすると
Pi to 10 decimal places: 3.1415926535
と表示され、
compute-pi 1000
とすると
Pi to 1000 decimal places: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
と表示される。クレートを読み込んで Rust プログラムの中から関数を使うこともできる。ドキュメント参照。
計算時間
MacBook Air (Apple M1, 16 GB)での計算時間を time コマンドで計測した。
- 1万桁: 15.85 ミリ秒
- 10万桁: 0.12 秒
- 100万桁: 1.51 秒
- 1000万桁: 25.30 秒
- 1億桁: 403.05 秒 (6.7 分)
- 2億桁: 14.67分
- 3億桁: 22.65 分
- 3億2000万桁: 23.52 分
- 3億3000万桁: 10時間で計算停止せず
Mac mini (Apple M1, 16 GB) ではこのようになった。
- 100万桁: 1.76 秒
- 1000万桁: 34.2 秒
- 1億桁: 498.80 秒 (8.3 分)
- 3億2000万桁: 28.48 分
- 3億3000万桁: 10時間で計算停止せず
このように、いずれの場合でも3億2000万桁までは30分以内に計算ができるが、3億3000万桁の計算はできなかった。メモリ割り当てに失敗してパニックで終了するのではなく、いつまでも計算が終了しないことから、メモリのスワッピングに時間がかかっているものと推察される。実際に、vm_stat で確認すると、Pageins と Pageouts が計算中に増えていく。16GBのメモリのマシン2台で同じ結果となったことから、16GBのメモリの場合はこれが限界の桁数であると考えられる。
なお、rug::float がu32で定義されていることから、メモリとは無関係に1,292,913,983桁が上限となっている。
他のプログラムとの比較
円周率計算の世界記録更新によく使われている y-cruncher で、同様に3億2000万桁を計算した。Mac mini で Docker の Debian コンテナから y-cruncher を起動して計算したところ、87秒で計算できた。MacBook では、なぜか2.2時間もかかった。このプログラムと y-cruncher の計算結果を比較したところ、3億2000万桁すべての数字が一致することが確認された。
Wikipedia のガウス・ルジャンドルのアルゴリズムの記事に掲載されている Python のプログラムをもとに、Paiza でブラウザ上で円周率の指定桁数を計算するプログラムを作成した。「入力」に桁数を入れて実行すると、指定桁数まで計算する（最終桁は丸める）。このプログラムで1万桁の計算はできたが、2万桁の計算は一定時間内に計算が終わらないため Timeout となった。Mac mini でこのプログラムを実行したところ、10万桁では11秒、100万桁では3分かかった。すなわち、Rust の compute-pi では、ガウス・ルジャンドルのアルゴリズムによる円周率100万桁の計算を Python プログラムの100倍の速度で実行できたということになる。