モンテカルロ法による、円周率の計算（C++）

　別記事〔下手な鉄砲の利用法〕で、BASICによるモンテカルロ法の計算例を紹介しました。しかし、計算回数が1000万回では、期待したほどの精度になっていないことがわかります。
　同じアルゴリズムを、C++を使ってプログラムに書いてみました。BASICに比べて、計算速度が速いと期待していました。でも、さすがに4億回の計算を10セットやらせたら、それなりの時間が掛かりましたね…＾＾；

　40億回計算し、(x^2+y^2)^(1/2)が1.0より小さかったのは、3,141,596,136回となっていますので、円周率のpi=3.141592654に近づいたことがわかります。

　でも、こんなに時間が掛かった割に、5桁の精度しか出ていません。モンテカルロ法というのは、なかなか難しいものだと思います…＾＾；

（by　心如）

─────

─────
// モンテカルロ法による円周率の計算
// ── Random（乱数）の応用例 ──
// 2012.04.08 ───── coded by 心如

#include <iostream>
#include <ctime>
using namespace std;

// 乱数の生成
class Random{
    int seed1, seed2, seed3;
    void SetSeed(int);
public:
    Random(){ SetSeed(1); }
    Random(int seed){ SetSeed( seed); }
    float rand();
    int rand(int n){ return (int)(rand() * n); }
};

// 乱数─────初期値の設定
void Random::SetSeed(int seed){
    seed1 = seed % 30269 + 1;
    seed2 = seed % 30307 + 1;
    seed3 = seed % 30323 + 1;
}

// 実数（0以上1.0未満）の一様乱数の発生
float Random::rand(){
    float fv;
    seed1 = (seed1 % 177) * 171 - (seed1 / 177) * 2;
    if(seed1 < 0) seed1 += 30269;
    seed2 = (seed2 % 176) * 172 - (seed2 / 176) * 35;
    if(seed2 < 0) seed2 += 30307;
    seed3 = (seed3 % 178) * 170 - (seed3 / 178) * 63;
    if(seed3 < 0) seed3 += 30323;
    fv = seed1 / 30269.0 + seed2 / 30307.0 + seed3 / 30323.0;
    while( fv >= 1.0) fv -= 1.0;
    return fv;
}

const Mn = 10000 * 10000;

// メイン・ルーチン
int main(){
    float x, y, r;
    Random ra((int)time(NULL));
   
    unsigned long s = 0;
    for(int j = 0;j < 10;j++){
        int cnt = 0;
        for(int i = 0;i < Mn * 4; i++){
            x = ra.rand();
            y = ra.rand();
            r = sqrt(x * x + y * y);
            if(r < 1.0) cnt++;
        }
        printf(" %d 回目: %d (%5.3f 秒)",
            j, cnt, (float )clock() / CLOCKS_PER_SEC);
        cout << "--" << ra.rand(100) << "--" << ra.rand() << endl;
        s += cnt;
    }
    cout << "s = " << s << endl;
}
─────

タグ：円周率 乱数