Mae向きなブログ

Mae向きな日記のブログ版。ようやくこちらに移行してきました。

モンテカルロ法

モンテカルロ法を使って円周率を求めるアルゴリズムは有名だけど,今まで実際にプログラムを作ってπを求めたことはなかった。
ちょっとやってみたところ,ちゃんと3.14に収束している。やっぱり,実際に動かしてみることは大切だと思う。

以下は,実行結果。sampleが点を打った数です。

sample =         10, pi = 2.8000
sample =        100, pi = 3.2000
sample =       1000, pi = 3.1440
sample =      10000, pi = 3.1676
sample =     100000, pi = 3.1453
sample =    1000000, pi = 3.1406
sample =   10000000, pi = 3.1419

以下は,サンプルプログラムです。

#include
#include
#include

#define SAMPLE 10000000

main()
{
  int cnt, max;
  int i;
  double x, y, pi;
  
  srand((unsigned)time(NULL));
  for(max = 10 ; max <= SAMPLE ; max *= 10){
    for (cnt = 0, i = 0;i < max;i++){
      x = (double)rand() / RAND_MAX;
      y = (double)rand() / RAND_MAX;
      if (x * x + y * y <= 1.0) cnt++;
    }
    pi = 4.0 * cnt / max;
    printf("sample = %10d, pi = %.4f\n",max,pi);
  }
}