こんにちは！今回は、Rustを使って円周率(π)を高速に計算する方法について紹介します。一般的なパソコンでも、わずか3秒程度で100万桁のπを計算できる「piday25」というプロジェクトを例に、その仕組みと実装方法を解説します。
はじめに
πの計算は、コンピュータの計算能力を評価するベンチマークとしても長く用いられてきました。現在では、様々なアルゴリズムによって数兆桁単位の計算も可能になっていますが、そのなかでもChudnovskyアルゴリズムは特に高速で効率的な方法として知られています。
今回紹介するpiday25は、このChudnovskyアルゴリズムにバイナリ分割法と並列処理を組み合わせることで、通常のデスクトップPCでも驚くべき速度でπを計算できる実装です。Intel Core i7（第10世代）のPCで、100万桁のπをわずか3秒程度で計算できるとのことです。
Chudnovskyアルゴリズムとは
まず、Chudnovskyアルゴリズムについて簡単に説明します。このアルゴリズムは1988年にChudnovsky兄弟によって発表されたもので、π計算の世界記録を達成するためにも使用されています。
数学的には以下の式で表されます：
1/π = (426880 * √10005) * Σ ((-1)^k * (6k)! * (13591409 + 545140134*k)) / ((3k)! * (k!)^3 * 640320^(3k))
この式は非常に速く収束するため、比較的少ない項数でも高精度の計算が可能です。たとえば、わずか14項ほどで100万桁近くの精度を得ることができます。
バイナリ分割法
Chudnovskyアルゴリズムをさらに高速化するために使われるのが「バイナリ分割法」です。これは総和計算など項が独立している場合に効率よく計算できる手法で、再帰的に計算を分割していきます。
バイナリ分割法の特徴：
- 計算を二分木のように分割して再帰的に処理
- 各項の計算が独立しているため並列処理に適している
- 計算量をO(n log n)に削減できる
piday25の実装解説
では、実際にpiday25の実装を見ていきましょう。コードは主に以下の部分から構成されています：
- コマンドライン引数の処理（main.rs）
- Chudnovskyアルゴリズムの実装（lib.rs）
- バイナリ分割と並列処理の仕組み
コアとなる処理の流れ
Rustコードの核心部分を簡略化すると、以下のような流れになります：
- binary_split関数で計算を再帰的に分割
- 複数のスレッドで並列に計算
- 計算結果をsum_up_par関数で集約
- 最終的にπの値を導出
バイナリ分割の実装
// 再帰的関数で、自分自身を2回呼び出して条件を満たすまで処理
fn binary_split(prec: &[i64; 6], j: i64, i: i64) -> (IBig, IBig, IBig) {
if i == j + 1 {
// 基底ケース：単一項の計算
let pab = IBig::from(-(6 * j - 5) * (2 * j - 1) * (6 * j - 1));
let qab = IBig::from(j).pow(3).mul(prec[3]);
let rab = IBig::from(prec[4] * j + prec[5]).mul(&pab);
(pab, qab, rab)
} else {
// 分割統治：範囲を半分に分けて再帰的に計算
let m = (j + i) / 2;
let (pam, qam, ram) = binary_split(prec, j, m);
let (pmb, qmb, rmb) = binary_split(prec, m, i);
// 結果を結合
let pab = pmb.mul(&pam);
let qab = qam.mul(&qmb);
let rab = ram.mul(qmb) + pam * rmb;
(pab, qab, rab)
}
}
この関数が再帰的に呼び出され、計算範囲を細かく分割していきます。
並列処理の実装
piday25では、crossbeam-channelクレートを使ってスレッド間の通信を行い、複数のCPUコアで並列に計算を進めています。
// 並列計算の部分（簡略化）
let m = (iterators as i64 + 1) / threads as i64;
let (s0, r) = unbounded();
// 複数のスレッドを起動して計算を並列化
for i in 1..threads {
let s = s0.clone();
thread::spawn(move || {
s.send((binary_split(prec, m * (i - 1) as i64, m * i as i64), i))
.unwrap();
});
}
各スレッドは割り当てられた範囲でバイナリ分割を実行し、結果をチャネルを通じてメインスレッドに送信します。
並列処理のポイント
このプログラムの並列処理の特徴として、以下の点が重要です：
- スレッド数は2のべき乗に制限：これは効率的な結果の集約に関係しています
- 各計算範囲の独立性：各スレッドは他に依存せず計算を完了できる
- 集約処理：sum_up_par関数で各スレッドの結果をツリー構造で効率的に集約
実行方法と性能
piday25を使用するには、まずRustのツールチェーンをインストールした上で、以下のコマンドを実行します：
git clone "https://github.com/elkasztano/piday25"
cd piday25
cargo build --release
コマンドライン引数
プログラムは次のようなコマンドライン引数をサポートしています：
-
-d, --digits：計算する桁数（デフォルト: 30,000） -
-t, --threads：使用するスレッド数（デフォルト: システムの並列処理能力） -
-v...：詳細度レベル -
-m, --measure-time：計算時間を表示 -
--silent：結果を出力しない（計算のみ） -
-o, --output：結果をファイルに保存
実行例
100万桁のπを計算し、最後の200桁だけを表示する例：
cargo run --release -- -d 1000200 -m -vvv | tail -c 200
この例では、計算時間も表示されます。
私のMacbook Pro (M2)では下記の結果でした
calculation of 1000200 digits of pi took 3.094 secs
プロジェクトの依存関係
piday25は以下の外部クレートに依存しています：
- dashu：任意精度の数値計算を行うためのライブラリ
- crossbeam-channel：並列計算用のチャネル実装
- clap：コマンドラインインターフェース
まとめ
piday25は、Chudnovskyアルゴリズムにバイナリ分割法が適用できることを見出し、そして並列処理を組み合わせることで、驚異的な速度でπを計算できるようにした素晴らしい実装です。
Rustでの並列処理の実装方法や、効率的な数値計算のテクニックを学ぶうえで良い教材でした。