モンテカルロ法で円周率を計算[C言語]

まぁやりつくされたネタではあるけれども，暇つぶしには丁度いいんじゃないかなと思って書くことにします。

1.モンテカルロ法とは

　簡単に言えば何かをシミュレーションをして，結果から求めたい値を探すというもの。1番わかりやすい例としてはサイコロ。サイコロが2個あります。2つとも同じ目が出る確率はいくつですか，という問題に対して実際に100回ぐらいサイコロふってみればわかるじゃないか，という話。そこから発展して，確率を求めることで円周率を計算したり数値的に積分したりできる。

2.具体的な計算法

　モンテカルロ法のWikipediaに良い図があるけど，一応一通り説明します。

　一辺が1の正方形の中に半径1の円が描かれています。ここで正方形の中に適当に砂を1粒落とすと，砂が円の内部に入る確率はいくつでしょう。

　正方形全体の面積は1，円の面積は半径1の円の4分の1の面積が描かれているからπ/4。したがって(確率) = (円の面積) ÷ (正方形全体の面積) だから答えはπ/4になる。

　じゃあ実際に適当に砂を，n粒落としたらm粒円の中に入ったとすると，

　

　と，こういう方法で近似的に円周率を求めることができる。

3.プログラム

　まぁ説明はあまり得意じゃないからね。プログラムで（お察しください）ということです。アルゴリズムは比較的単純。砂粒を落とすという表現を，一辺1の正方形の中にランダムに点をおくということに言い換える。

　m=0から始めて，

　「x,yに0から1の乱数を代入→円の中ならmを1増やす」という作業をn回繰り返す。

　円の中なら原点からの距離が1以下(未満)であれば良い。

/**
 * main.c
 * 2015/09/29
 * モンテカルロ法で円周率を求めるプログラム
 */

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

/* 試行回数 */
#define N 100000

/* ランダムに点を決め，円の中なら1，そうでなければ0を返す */
int step(){
   double x,y;
   
   /* x,yに0から1の乱数を代入する */
   x = (double)rand() / ( RAND_MAX + 1);
   y = (double)rand() / ( RAND_MAX + 1);
   
   /* 原点からの距離の二乗が1未満なら円の中 */
   if( x * x + y * y < 1 )
      return 1;
   
   return 0;
}

int main(void){
   int i;
   int m = 0;
   
   /* 乱数の種をまく */
   srand(time(NULL));
   
   for( i = 0 ; i < N ; i++ ){
      /* 点が円の中なら */
      if( step() == 1 ){
         /* mを1増やす */
         m++;
      }
   }
   
   /* 結果を出力する */
   printf("pi = %lf\n", 4.0 * m / N );
      
   return 0;
}

出力結果はpi = 3.144840と，まぁそこそこあってる。でも物足りない。計算した気がしない。円周率をゴリゴリ計算するプログラムはこっち。まぁこれも遅いけど一応任意の桁数を計算できます。

追記：Processingで可視化した

だいたいみんな書いたことがあるモンテカルロ法で円周率を求めるプログラムを可視化してみた…… pic.twitter.com/Znnuxh34HX

— もりとにー (@TonyMooori) 2016, 2月 4