SSブログ

モンテカルロ法による、円周率の計算(高速化) [C++]

 別記事〔モンテカルロ法による、円周率の計算(C++)〕で、C++で描いたプログラムを紹介しました。しかし、コンパイラ言語であるC++を使っても、40億回も計算を繰り返すと、かなりの時間がかかることがわかりました。
 試しに、乱数発生部分をもっと簡単なものに書き換えてみましたので紹介します。別記事で、1,314秒も掛かっていた計算が、約291秒で出来ました。

 40億回計算し、r=sqr(x^2+y^2)が1.0未満になったのは、3,141,718,305回という結果になりました。しかし、残念ですが、簡易タイプの乱数では、計算速度は速くなっていますが、計算精度は低下しました。

 今回は、高速であっても一様な乱数を生成できないと、計算精度は低下するという、悪いプログラムの見本となってしまいました…^^;

(by 心如)

─────
CPP cd010.JPG


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

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

// 乱数(0~32767の整数)
class Random{
    unsigned long nxt;
public:
    Random(){ nxt = 1; }
    Random(int sv){ nxt = (unsigned long) sv; }
    int i(){
        nxt = nxt * 1103515245 + 12345;
        return (int)((nxt / 65536) % 32768);
    }
    float f(){ return (float)i() / 32768.0; }
};

const Mn = 100000000;

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

【追記】
 ちなみに、同じ計算を、C++の組み込み関数rand()を使ってやってみました。計算速度はほぼ同じですが、精度はほんの少し悪いか同じという結果になっていると思います。
CPP cd011.JPG

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

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

// 乱数(0~32767の整数)
class Random{
public:
    Random(){ srand(time(NULL)); }
    int i(){ return rand(); }
    float f(){ return (float)i() / 32768.0; }
};

const Mn = 100000000;

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


タグ:円周率 乱数
nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。