C 言語でモンテカルロ法を用いた円周率 (π) の近似計算方法

はじめに

この実験では、C プログラミング言語を用いてモンテカルロ法を用いて円周率 π の近似値を求める方法を学びます。まず、単位正方形内にランダムな点を生成し、そのうち四分円内に含まれる点の数をカウントすることで、π の近似値を計算します。最後に、計算された近似値を出力します。

この実験では、π の推定のためのモンテカルロ法の主要なステップ、すなわちランダムな点の生成、四分円内にある点のカウント、そして総点数のうち四分円内にある点の比率を用いて近似値を計算する手順を網羅しています。この実験は、C プログラミング言語を用いて微積分や解析幾何学の概念を実用的に適用する例を示しています。

単位正方形内のランダムな点の生成

このステップでは、C プログラミング言語を用いて単位正方形内にランダムな点を生成する方法を学びます。これは、円周率 π の近似値を求めるモンテカルロ法において重要な部分です。

まず、~/project ディレクトリに新しい C ファイルを作成し、ランダムな点の生成を実装しましょう。

cd ~/project
nano random_points.c

次に、ランダムな点の生成を行うコードを記述します。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define NUM_POINTS 10000

int main() {
    // 乱数生成器のシードを設定
    srand(time(NULL));

    // 単位正方形内のランダムな点を生成
    for (int i = 0; i < NUM_POINTS; i++) {
        double x = (double)rand() / RAND_MAX;
        double y = (double)rand() / RAND_MAX;

        printf("点 %d: (%.4f, %.4f)\n", i + 1, x, y);
    }

    return 0;
}

コードの詳細を説明します。

srand(time(NULL)) は、現在の時刻で乱数生成器のシードを設定します。
(double)rand() / RAND_MAX は、0 から 1 の間のランダムな数を生成します。
10,000 個のランダムな点を単位正方形 (0,0) から (1,1) に生成します。

プログラムをコンパイルして実行します。

gcc random_points.c -o random_points
./random_points

実行結果の例：

点 1: (0.7234, 0.5678)
点 2: (0.2345, 0.9876)
...
点 10000: (0.1122, 0.3344)

四分円内点のカウントと円周率近似

このステップでは、前のプログラムを修正して四分円内にある点の数をカウントし、その情報を使って円周率 π の近似値を求めます。

random_points.c ファイルを更新しましょう。

cd ~/project
nano random_points.c

以前のコードを以下の実装に置き換えます。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define NUM_POINTS 100000

int main() {
    // 乱数生成器のシードを設定
    srand(time(NULL));

    int points_inside_circle = 0;

    // ランダムな点を生成し、四分円内にある点のカウント
    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) {
            points_inside_circle++;
        }
    }

    // 円周率の近似値を計算
    double pi_approximation = 4.0 * points_inside_circle / NUM_POINTS;

    printf("総点の数：%d\n", NUM_POINTS);
    printf("四分円内にある点の数：%d\n", points_inside_circle);
    printf("円周率の近似値：%.6f\n", pi_approximation);
    printf("実際の円周率の値：   %.6f\n", M_PI);
    printf("誤差：       %.6f\n", fabs(pi_approximation - M_PI));

    return 0;
}

数学ライブラリを使用してプログラムをコンパイルします。

gcc random_points.c -o random_points -lm
./random_points

実行結果の例：

総点の数: 100000
四分円内にある点の数: 78540
円周率の近似値: 3.141600
実際の円周率の値:    3.141593
誤差:        0.000007

変更点の詳細を説明します。

精度向上のため、点の数を増やしました。
四分円内にある点のカウントロジックを追加しました。
式：π ≈ 4 * (四分円内にある点の数) / (総点の数) を使用しました。
実際の円周率の値との比較を含めました。

近似値の表示

この最終ステップでは、円周率の近似値を求めるプログラムを拡張し、結果を表示する関数を作成して出力形式を改善します。

random_points.c ファイルを修正しましょう。

cd ~/project
nano random_points.c

新しい表示関数と改良された出力でコードを更新します。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define NUM_POINTS 1000000

// 近似結果を表示する関数
void print_pi_approximation(int total_points, int points_inside_circle) {
    double pi_approximation = 4.0 * points_inside_circle / total_points;

    printf("円周率近似結果\n");
    printf("=====================\n");
    printf("生成した総点の数：       %d\n", total_points);
    printf("四分円内にある点の数： %d\n", points_inside_circle);
    printf("近似された円周率の値：         %.8f\n", pi_approximation);
    printf("実際の円周率の値：               %.8f\n", M_PI);
    printf("絶対誤差：          %.8f\n", fabs(pi_approximation - M_PI));
    printf("近似精度：       %.4f%%\n",
           (1 - fabs(pi_approximation - M_PI) / M_PI) * 100);
}

int main() {
    // 乱数生成器のシードを設定
    srand(time(NULL));

    int points_inside_circle = 0;

    // ランダムな点を生成し、四分円内にある点のカウント
    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) {
            points_inside_circle++;
        }
    }

    // 近似結果を表示
    print_pi_approximation(NUM_POINTS, points_inside_circle);

    return 0;
}

プログラムをコンパイルして実行します。

gcc random_points.c -o random_points -lm
./random_points

実行結果の例：

円周率近似結果
=====================
生成した総点の数:        1000000
四分円内にある点の数:  785398
近似された円周率の値:          3.14159200
実際の円周率の値:                3.14159265
絶対誤差:           0.00000065
近似精度:        99.9998%

主な改善点：

専用関数 print_pi_approximation() を作成しました。
より高い精度を得るために、点を 1,000,000 個に増やしました。
出力形式をより詳細に改善しました。
近似精度のパーセンテージを含めました。

まとめ

この実験では、最初に C プログラミング言語を用いて単位正方形内にランダムな点を生成する方法を学びました。これは、円周率 π の近似値を求めるモンテカルロ法の重要な部分です。乱数生成器のシードを現在時刻で初期化し、rand() 関数を使って単位正方形 (0,0) から (1,1) の範囲に 10,000 個のランダムな点を生成しました。

次に、プログラムを修正して四分円内にある点の数をカウントし、その情報を使って円周率 π の近似値を求めました。100,000 個のランダムな点を生成し、それらが四分円内に含まれるかどうかをチェックしました。円内にある点の数と総点数の比率を計算することで、円周率 π の値を推定することができました。