Seki's note

Rust による円周率計算プログラム

著者：関 勝寿
公開日：2024年3月17日
キーワード： 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万桁を計算した。以下に、いくつかの参考サイトを示す。

使い方

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 コマンドで計測した。

Mac mini (Apple M1, 16 GB) ではこのようになった。

このように、いずれの場合でも3億2000万桁までは30分以内に計算ができるが、3億3000万桁の計算はできなかった。メモリ割り当てに失敗してパニックで終了するのではなく、いつまでも計算が終了しないことから、メモリのスワッピングに時間がかかっているものと推察される。実際に、vm_stat で確認すると、Pageins と Pageouts が計算中に増えていく。16GBのメモリのマシン2台で同じ結果となったことから、16GBのメモリの場合はこれが限界の桁数であると考えられる。

他のプログラムとの比較

円周率計算の世界記録更新によく使われている 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倍の速度で実行できたということになる。