C言語でモンテカルロ法を用いた円周率の計算
はじめに
円周率（π）は数学や物理学、工学で重要な定数です。πの計算方法にはさまざまな手法がありますが、本記事では モンテカルロ法 を使って円周率を求める方法をC言語で実装します。
モンテカルロ法とは？
モンテカルロ法は、乱数を用いた確率的手法の一つで、シミュレーションを通じて数値を近似的に求めることができます。円周率の計算においては、ランダムな点を使って 円の面積と正方形の面積の比率 を利用してπを推定します。
1. モンテカルロ法による円周率の求め方
モンテカルロ法を使って円周率を求める基本的なアイデアは以下の通りです。
- 1×1 の正方形の中に ランダムな点（x, y） を大量に生成する
- それぞれの点が、原点中心の 半径1の円の内部 に入っているかを判定する（円の方程式：x² + y² ≤ 1）
- 円の内部に入った点の割合を計算し、円周率を近似する
円の面積は π × r²（ここでは r = 1 なので π × 1² = π）
正方形の面積は 1 × 1 = 1 なので、
円の内部に入る確率は π/4 となります。
よって、以下の式で円周率を求めます。
[ \pi ≈ 4 × \frac{\text{円の内部に入った点の数}}{\text{全体の点の数}} ]
2. C言語での実装
必要なライブラリ
stdlib.h（乱数生成）stdio.h（入出力）math.h（平方根計算）
コードの実装
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> // モンテカルロ法による円周率の計算 double monte_carlo_pi(int num_points) { int count_inside_circle = 0; for (int i = 0; i < num_points; i++) { // 乱数生成 (0 ～ 1 の範囲) double x = (double)rand() / RAND_MAX; double y = (double)rand() / RAND_MAX; // 円の内側にあるか判定 if (x * x + y * y <= 1.0) { count_inside_circle++; } } // 円周率の推定値を計算 return 4.0 * count_inside_circle / num_points; } int main() { // 乱数の種を設定（現在時刻を利用） srand(time(NULL)); // 試行回数を指定 int num_points = 1000000; // 100万点 // πを計算 double estimated_pi = monte_carlo_pi(num_points); // 結果を表示 printf("モンテカルロ法による円周率の推定値: %.6f\n", estimated_pi); return 0; }
3. コードの解説
1. 乱数の生成
double x = (double)rand() / RAND_MAX; double y = (double)rand() / RAND_MAX;
rand() は 0 から RAND_MAX までの整数を返します。
これを RAND_MAX で割ることで 0.0 から 1.0 までの実数を取得できます。
2. 円の内部に入っているか判定
if (x * x + y * y <= 1.0) { count_inside_circle++; }
原点 (0,0) からの距離が 1.0 以下であれば、円の内部にあると判定します。
3. πの近似計算
return 4.0 * count_inside_circle / num_points;
円の内部に入った点の割合に 4.0 を掛けることでπの推定値を求めます。
4. 乱数の種の設定
srand(time(NULL));
乱数の種を time(NULL) で設定することで、プログラム実行ごとに異なる乱数を生成します。
4. 実行結果
試行回数を増やすと、より正確な円周率が求められます。
例えば、num_points = 1,000,000（100万回）の場合、実行すると以下のような結果が得られます。
モンテカルロ法による円周率の推定値: 3.141593
試行回数を変えた場合の精度の違い： | 試行回数 | 近似値（例） | |---------|------------| | 1,000 | 3.148000 | | 10,000 | 3.141200 | | 100,000 | 3.141540 | | 1,000,000 | 3.141593 |
試行回数が多いほど精度が向上することが分かります。
5. さらなる改良
並列化による高速化
OpenMP を使用してループを並列処理することで、高速化できます。
#include <omp.h> double monte_carlo_pi_parallel(int num_points) { int count_inside_circle = 0; #pragma omp parallel for reduction(+:count_inside_circle) for (int i = 0; i < num_points; i++) { double x = (double)rand() / RAND_MAX; double y = (double)rand() / RAND_MAX; if (x * x + y * y <= 1.0) { count_inside_circle++; } } return 4.0 * count_inside_circle / num_points; }
このコードでは、#pragma omp parallel for reduction(+:count_inside_circle) を使って、各スレッドが並行して計算を行い、結果を合計するようにしています。
6. まとめ
本記事では、C言語を用いて モンテカルロ法による円周率の計算 を実装しました。
ポイント
- モンテカルロ法を使って乱数を用いた円周率の推定を行った
- 乱数を生成し、円の内部に入る割合を使ってπを求めた
- 試行回数を増やすと精度が向上することを確認した
- 並列処理を導入すれば高速化できる
