Rustで円周率計算&収束の速さを検証してみた（単一CPU）

せっかくRaspberry Piクラスタを構築したので、並列計算を試してみたい、ということで円周率を計算してみることにしました。ただC++でプログラムを書いても面白くないので、以前から関心があったRustの学習を兼ねて実装してみました。

Rust環境構築

Rustのインストール

curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh

もしくはHomebrewの環境であれば

brew install rust

でインストールできます。

依存パッケージの追加

円周率を大きな桁数で計算するためには、任意精度で計算を行える必要があります。今回はrugを用いて任意精度計算を実現します。rugは、GMP, MPFR, MPCをRustから使えるようにしたものです。整数、有理数、浮動小数、複素数の演算を任意精度で行うことができます。以下のコマンドを実行するとrugを追加できます。

cargo add rug

poetryと同じような使い方かな？という印象です。

コードのフォーマット

プログラムを書くと自動で整形したくなるもので、フォーマッタについても調査しました。どうやらフォーマッタも付属のようで、以下のコマンドでフォーマットできました。便利〜。

cargo fmt

円周率の計算

円周率の計算方法については、こちらの記事を参考にしました。今回は64ビット浮動小数点型での円周率との差を基準に、以下の計算方法を用いて円周率を求める場合の収束の速さを比較します。

実装

参考ウェブサイトの数式をただプログラムに直しただけなので、特殊な工夫等はしていません。ソースコードは並列版を作成してからまとめてgithubで公開しようかな…。

実験結果

以下は縦軸が2乗誤差、横軸が何項目までの和を取るかを表示した両対数グラフです。収束が遅いとされているライプニッツの式やWallis級数では1,000,000項目まで計算しても2乗誤差が$10^{-5}$程度しか下がらないのに対し、幾何平均ベースでは横軸の値が小さくてもかなり2乗誤差が小さいことがわかります。想像していた以上に収束が遅いことがわかりました。

並列化に向けて

今回実装したsequentialなアルゴリズムはそのまま並列化することができません。並列実行するためには、Binary Splitting methodやDRM法を用いる必要があります。どうせなら速度にも拘りたいので、現在最も速いとされているChudnovskyの公式を用いた並列化プログラムを組みたいと考えています。

参考

• 円周率を求める式
• 128ビット符号付き整数の最大値は素数 – Rustで任意精度整数演算
• 多倍長計算を行う (rug crate)
• 円周率.jp
• 100日後にRustをちょっと知ってる人になる: [Day 69]コード整形: rustfmt