🦀

2600億回のモンテカルロを14秒で：Rust×rayon×SIMDで並列化してみた

に公開
Rust
performance
SIMD
multithreading
zennfes2025free
tech

1. 導入 🦀

モンテカルロ法はシンプルながらも膨大な試行回数を必要とするアルゴリズムです。数十億、数百億といった規模になると、素朴な実装では到底現実的な時間で終わりません。
今回はこのモンテカルロ法で円周率を推定しようと試みました。
今回取り組んだ末に得たものは「2600億回の試行を14秒で終わらせる」でした。Rustの安全性と低レベル最適化の両立、そして並列化のためのrayon、さらにCPUの演算能力を引き出すSIMDを組み合わせることで、CPUフレンドリーかつの余力を使い切る実装を目指しました。

2. モンテカルロ法による円周率推定

モンテカルロ法で円周率を推定するためには、以下の方法を使います。

  1. 中心が(0, 0)の半径1の円の第一象限を用意する。
  2. 0 \le x \le 1, 0 \le y \le 1でランダムに点(x, y)を取る。
  3. その点が円の中に入っているか、すなわちx^2 + y^2 \le 1を判定する。
  4. 2,3を繰り返す。
  5. 円の中に入った点の割合をpとすると、円周率の近似値は
    \pi \approx 4p = 4 \times \frac{\text{円の中に入った点の数}}{\text{総試行回数}}

モンテカルロ法は精度が上がるのが非常に遅く、実行回数をNとすると、信頼できる桁数は\frac{1}{\sqrt{N}}\text{桁}となるそうです。

3. 設計方針

今回は「まず素朴に書いてから改善する」のではなく、最初から最適化を前提に設計しました。
設計の基本方針は以下の3点です：

  • CPUコアを余すことなく使う(マルチスレッド化)
  • 各コアの演算器をフル稼働させる(SIMD化)
  • CPUにとって実行しやすい形にする(最適化)

4. 並列化の実装(rayon)

Rustで並列化を行う際、rayonは非常に強力です。イテレータをpar_iterに置き換えるだけで、複数スレッドに処理を分散できます。
ただし「ただ並列にすれば速い」というわけではありません。各スレッドが同じデータを操作するためのAtomic操作と呼ばれるものを使うと、CPUが同期されて遅くなるため、回数を極力下げて運用しなければなりません。
また、すべてのスレッドを使用して実行すると、OSや他プロセスとのバッティングが起き、実行速度が低下するため、1スレッドだけ残しておくのが賢明かもしれません。

5. 乱数生成器の最適化

このモンテカルロ法を高速で処理するにあたって、最もボトルネックな要素は乱数の生成です。randクレートが提供するthread_rngによる生成は高精度な代わりに致命的に実行速度が遅いです。そのため、今回の用途には適しているとは言えません。今回はXorShiftというXorとShiftのみで構成されている乱数生成器を使用します。これは乱数としての精度が低い代わりに、非常に高速で乱数を作ることが可能です。また、のちに紹介するSIMDとの相性も抜群です。

6. SIMD最適化の実装

SIMD(Single Instruction, Multiple Data)は、同じ演算を複数データに同時適用できる仕組みです。Rustでは将来的に、std::simdを利用することで比較的直感的にベクトル化が可能になる予定です。しかし、今回はstd::arch::x86_64の命令を使いました。これらはCPUの命令セットレベルの低レベルな命令で、unsafeを要求します。
今回の実装では、乱数生成やループ内の演算処理をSIMD化し、単精度浮動小数の利用により、1命令で8回の試行を同時に進められるようにしました。これにより、単純なループ実装に比べて演算密度が大幅に向上します。

7. 命令の最適化

乱数のスケーリング

プログラムを高速化するためには、まず、余分な命令を減らすことが重要です。今回の例では[0,1)の乱数を生成する必要がありますが、単純に乱数を整数で取り出してから「最大値で割る」といったスケーリングを行うと、除算命令がボトルネックになります。除算はCPUにとって最も遅い命令のひとつだからです。
そこで注目するのがIEEE754単精度浮動小数点数の内部表現です。単精度浮動小数点は以下のように構成されています。

  • 符号部：1ビット
  • 指数部：8ビット
  • 仮数部：23ビット

ここで「指数部に127を設定すると、その数は必ず[1,2)の範囲を表す」という性質を利用します。つまり、乱数ビットを仮数部に直接埋め込み、指数部を127に固定することで、[1,2)の乱数を一発で生成できるのです。
さらに、この値から単純に1を引けば[0,1)の乱数が得られます。必要な命令は「指数部をセットする処理」と「引き算」のわずか2つだけとなり、割り算を使ったスケーリングよりもはるかに高速に乱数を生成できます。
ちなみに、今回の方法では、浮動小数点の符号ビットもランダムに決まります。
ただし、モンテカルロ法で円の内部に入っているかどうかを判定するときは
x^2 + y^2 \leq 1
という形で二乗を使います。
そのため、符号が正でも負でも結果は同じになり、符号ビットは完全に無視できるのです。

円内の判定

円の内外判定は
x^2 + y^2 \leq 1
ただし、浮動小数点の比較命令はNaNの扱いや例外パスを考慮する必要があり、フラグ依存や分岐予測ミスによって実効スループットが低下する場合があります。そこで、比較命令を避けて減算＋符号ビット判定に置き換えることで、分岐を生まずに高速に処理できるようにしました。
式を変形すると：
x^2 + y^2 - 1 \leq 0
つまり「負かどうか」を判定すればよいということになります。浮動小数点数の先頭1ビットは符号部なので、単に符号ビットを見るだけで判定可能です。

  • 和から1を引く
  • 先頭ビットが1なら負
  • これをビット列として回収し、_popcnt()で一気に数える

こうすることで、比較命令を使わずにSIMD + ビット演算で高速に円内判定ができます。

8. 乱数の再利用

本来であれば、モンテカルロ法では毎回新しい乱数を使うことで独立性を保つのが理想です。しかし、今回の実装において最も重い処理は乱数生成でした。もし乱数を再利用できれば、その分だけ実行回数を増やせるため、時間当たりの試行回数がほぼ倍増します。
今回検討したのは、x座標に使った乱数をy座標にも再利用するという方法です。

SIMDと乱数の並び替え

今回の実装では SIMD によって同時に8つの乱数を生成しています。つまり常に
a, b, c, d, e, f, g, h
という8つの乱数が存在しています。

  • まずこれらをx座標として利用する
  • その後、並び替えてy座標として再利用する

こうすれば、乱数生成の回数は半分で済みます。

独立性の低下とその影響

もちろん、この方法では乱数の独立性が下がります。例えば：

  • 点A = (a, b)
  • 点B = (b, c)
  • 点H = (h, a)

このように定義すると、bが変化すると点Aのy座標と点Bのx座標が同時に変化します。
ただし、すべての点について、a や c など他の乱数が絡むため、完全に同じ挙動を取ることはありません。

トレードオフ

  • デメリット: 乱数の独立性が下がる
  • メリット: 乱数生成コストが半減し、実行回数が大幅に増える

モンテカルロ法では「実行回数こそが精度を決める」ため、この程度のランダム性低下は必要経費と割り切ることもできます。

9. rayon × SIMD の統合設計

rayonによるスレッド並列と、SIMDによるデータ並列を組み合わせることで、CPUを「横方向(スレッド)」と「縦方向(SIMDレーン)」の両面から使い切ることができます。

10. コンパイラフラグ

これがめちゃくちゃ重要です。きちんと設定したかどうかで、速度が二倍以上変わってくる場合もあります。私は以下のように設定しました。

# Cargo.toml
opt-level = 3        # O3相当の最適化
lto = true           # リンカ時最適化（Link Time Optimization）
codegen-units = 1    # 1スレッドでコード生成 → 最適化が効きやすい
panic = "abort"      # panic時にスタック展開せず即終了
strip = true         # デバッグ情報を削除してバイナリを軽量化

# .cargo/config.toml
rustflags = ["-C", "target-cpu=native"] # 実行環境のCPUに最適化された命令を使う

11. 実行結果

これからベンチマークの実行結果を紹介します。ただし、環境が異なると数値も変わってくるため、まずは私の実行環境を示しておきます。

  • CPU: Intel Core Ultra 5 125U(14スレッド)
  • RAM: DDR5 16GB
  • OS / Toolchain: Windows 11 + Rust 1.87.0

このCPUはモバイル向けの省電力モデルですが、AVX2命令セットを搭載しているため、256bit幅でのSIMD演算が可能です。ただしノートPCでの実行のため、数十秒の使用でサーマルスロットリングが発生し、常時フル稼働は難しい点に注意が必要です。

以下が結果になります。

interior:204203538220 total:260000000000 pi:3.1415928956923076 time: 14.5927911s

interior:204203778888 total:260000000000 pi:3.141596598276923 time: 14.8603206s

interior:204203707828 total:260000000000 pi:3.141595505046154 time: 14.1803344s

実行速度

  • 総試行回数: 2600億回
  • 実行時間: 約14.5秒
  • 1秒あたりの試行回数:
    \frac{260{,}000{,}000{,}000}{14.5} \approx 1.79 \times 10^{10} \ \text{回/秒}
  • → 毎秒約180億回の試行が可能でした。
  • 1スレッド当たりの処理量は約14億回となります。

精度

  • 推定された円周率はいずれも3.14159…まで安定して出力。
  • 理論的には、モンテカルロ法の誤差は\frac{1}{\sqrt{N}}で収束します。
  • 今回のN = 2.6 \times 10^{11}に対しては、
    \frac{1}{\sqrt{N}} \approx 2 \times 10^{-6}
    となり、実際に得られた誤差（10^{-6}オーダー）と一致しており、理想的な結果が得られました。
  • すなわち、乱数の再利用は適切な範囲でならば、結果に対してほとんど影響をもたらさないと言えるでしょう。

12. まとめと展望

今回の取り組みでは、私の環境に合わせた最適化を徹底することで、1秒あたり約180億回という膨大な試行を実現することができました。これは、モバイル向けCPUであっても工夫次第でここまで性能を引き出せる、という一つの証明になったと思います。
今後の展望としては、さらなる高速化の余地がいくつか考えられます。

  • AVX-512
    AVX2の2倍の幅を持つベクトル演算に加え、多様な命令セットが利用可能です。特に、固定小数点（Q2.14形式など）での計算に有効な「積の上位ビットを直接取得する命令」が存在するため、単純に2倍どころか、さらに効率的な実装が期待できます。
  • GPUによる並列化
    モンテカルロ法は「大量の独立した試行」を繰り返すアルゴリズムであり、GPUの得意分野そのものです。数千〜数万スレッドを同時に走らせることで、CPUでは到底届かないスループットを実現できるでしょう。

コード例とリポジトリ

今回紹介した実装の全コードは以下のリポジトリに公開しています。
興味のある方はぜひcloneして試してみてください。
実行回数の変更方法などについては、リポジトリのREADMEをご覧ください。

実行方法

git clone https://github.com/yua134/montecarlo-pi.git
cd montecarlo-pi
cargo build --release
cargo run --release

Discussion