モンテカルロ法を用いて円周率を確率的に求めるC言語のコード

 はい、モンテカルロ法を用いて円周率を確率的に求めるC言語のコードを提供します。

---

円周率を確率で求めるCコード

このCコードもPython版と同様に、モンテカルロ法を使って円周率を近似します。

基本的な考え方:

1. 一辺の長さが1の正方形の中に、原点を中心とする半径1の円の4分の1（扇形）を考えます。

2. 正方形の中にランダムに点を打ちます。

3. その点が扇形の中に入った回数を数えます。

4. 「扇形の中に入った点の数 / 全体の点の数」が、扇形の面積と正方形の面積の比 (pi/4) に近似することを利用して、円周率を推定します。

C

#include <stdio.h>   // printfなどの入出力関数用
#include <stdlib.h>  // rand, srandなどの乱数関数用
#include <time.h>    // time関数用（乱数のシード設定に利用）
#include <math.h>    // sqrt, powなどの数学関数用 (ここではpowのみ使用)

// モンテカルロ法を用いて円周率を推定する関数
double estimate_pi(long long num_darts) {
    long long inside_circle_count = 0; // 円の中に入った点の数
    long long i;
    double x, y; // 点の座標
    double distance_squared; // 原点からの距離の2乗

    // 乱数のシードを設定
    // これにより、プログラムを実行するたびに異なる乱数列が生成されます
    srand((unsigned int)time(NULL));

    for (i = 0; i < num_darts; i++) {
        // 0.0から1.0の範囲でランダムなx座標とy座標を生成
        // RAND_MAXはrand()が返す最大値
        x = (double)rand() / RAND_MAX;
        y = (double)rand() / RAND_MAX;

        // 原点からの距離の2乗を計算 (x^2 + y^2)
        distance_squared = pow(x, 2) + pow(y, 2);

        // 距離が半径1以内（つまり距離の2乗が1以内）であれば円の中に入ったと判定
        if (distance_squared <= 1.0) {
            inside_circle_count++;
        }
    }

    // 円周率の推定値を計算
    // (円の中に入った点の数 / 全体の点の数) * 4
    return 4.0 * ((double)inside_circle_count / num_darts);
}

int main() {
    // 試行回数を設定
    long long num_darts_1 = 1000;
    long long num_darts_2 = 100000;
    long long num_darts_3 = 1000000;
    long long num_darts_4 = 10000000; // さらに多くの矢を試す

    double pi_approx_1, pi_approx_2, pi_approx_3, pi_approx_4;

    printf("--- 円周率のモンテカルロ推定 ---\n");

    // 各試行回数での円周率の推定
    pi_approx_1 = estimate_pi(num_darts_1);
    printf("矢の数 %lld 本の場合の円周率の推定値: %.8f\n", num_darts_1, pi_approx_1);

    pi_approx_2 = estimate_pi(num_darts_2);
    printf("矢の数 %lld 本の場合の円周率の推定値: %.8f\n", num_darts_2, pi_approx_2);

    pi_approx_3 = estimate_pi(num_darts_3);
    printf("矢の数 %lld 本の場合の円周率の推定値: %.8f\n", num_darts_3, pi_approx_3);

    pi_approx_4 = estimate_pi(num_darts_4);
    printf("矢の数 %lld 本の場合の円周率の推定値: %.8f\n", num_darts_4, pi_approx_4);

    // 実際の円周率と比較
    printf("\n実際の円周率 (M_PI): %.8f\n", M_PI); // M_PIはmath.hで定義されている円周率

    return 0; // プログラムが正常終了したことを示す
}

---

コードの解説

1. ヘッダーファイルのインクルード:

   • stdio.h: printf() などの標準入出力関数を使うために必要です。

   • stdlib.h: rand() (乱数生成) や srand() (乱数シード設定) など、汎用的な関数を使うために必要です。

   • time.h: time() 関数を使うために必要です。これは srand() の乱数シードとして現在の時刻を使い、プログラムを実行するたびに異なる乱数列が生成されるようにします。

   • math.h: pow() (べき乗) 関数を使うために必要です。また、M_PI (円周率の定数) も通常このヘッダーで定義されています。

2. estimate_pi(long long num_darts) 関数:

   • long long num_darts: 投げる矢の総数です。大きな数になる可能性があるため long long 型を使用しています。

   • srand((unsigned int)time(NULL));: 乱数ジェネレータのシードを設定します。time(NULL) は現在の時刻を秒単位で返し、これをシードとすることで、プログラムを複数回実行しても異なる結果（ランダムな点の配置）が得られます。この行は、estimate_pi関数が複数回呼び出される場合でも、通常はmain関数内など、プログラムの開始時に一度だけ呼び出すのが一般的です。 今回は関数内で乱数生成を行うため、関数内に記述しています。

   • x = (double)rand() / RAND_MAX;: rand() 関数は0から RAND_MAX (システム定義の最大値) までの擬似乱数を整数で生成します。これを RAND_MAX で割ることで、0.0から1.0までの範囲の浮動小数点数に変換し、x座標として使用します。y も同様に生成します。

   • distance_squared = pow(x, 2) + pow(y, 2);: 生成された点 (x, y) の原点からの距離の2乗を計算します。pow(x, 2) は x2 を意味します。

   • if (distance_squared <= 1.0): 計算された距離の2乗が1.0以下であれば、その点が半径1の円（の4分の1）の内部または境界上にあると判断し、inside_circle_count をインクリメントします。

   • return 4.0 * ((double)inside_circle_count / num_darts);: 円の中に入った点の割合を計算し、それを4倍して円周率の推定値として返します。(double) でキャストすることで、整数同士の割り算ではなく浮動小数点数での割り算が行われるようにします。

3. main() 関数:

   • 様々な num_darts の値（試行回数）を設定し、それぞれのケースで estimate_pi 関数を呼び出して円周率を推定します。

   • printf() を使って結果を表示します。%.8f は、小数点以下8桁まで表示することを意味します。

   • M_PI は math.h で定義されている円周率の定数で、推定値と比較するために使用します。

---

コンパイルと実行方法

1. 上記のコードを monte_carlo_pi.c などのファイル名で保存します。

2. C言語コンパイラ（例えばGCC）を使ってコンパイルします。ターミナル（コマンドプロンプト）を開き、以下のコマンドを実行します。

   Bash

   gcc monte_carlo_pi.c -o monte_carlo_pi -lm
   

   • -o monte_carlo_pi: 出力ファイル名を monte_carlo_pi と指定します。

   • -lm: math.h 内の数学関数（pow() など）を使用する場合に必要なオプションです。

3. コンパイルが成功すると、実行可能ファイル (monte_carlo_pi または monte_carlo_pi.exe ) が生成されます。以下のコマンドで実行します。

   Bash

   ./monte_carlo_pi
   

矢の数を増やすほど、計算に時間がかかりますが、円周率の推定値が実際の値に近づくことが確認できるはずです。