適当にモンテカルロ法で円周率を計算してみた
はてなのシンタックスハイライト機能のテストもしたいなぁとか思ってモンテカルロ法による円周率計算を行うプログラムを書いてみた.
乱数発生
この手のシミュレーションには欠かせないのが乱数.
しかし, .NET標準の乱数発生器Randomクラスは
クヌース先生が作ったアルゴリズムを実装しているやつらしいけど,
生成される乱数の質がどうもよろしくないらしい.
(ここ参照: 良い乱数・悪い乱数)
こういう部分に無駄にこだわってしまう自分は
Randomクラスをそのまま使う気にはならないなぁ…
ということで質の良い乱数を生成してくれることで定評のある
メルセンヌツイスターを使ってみることにした.
## メルセンヌツイスターについて詳しく知りたい人へ
## いちいち想像するのも面倒だ, てめえで勝手に想像しろ(ベジータ風)
## 作者の一人である松本眞さんが
## 説明してくれているサイトがあるので是非そちらをどうぞ
## Mersenne Twister: A random number generator (since 1997/10)
メルセンヌツイスターの実装
面倒くさかったので
http://takel.jp/mt/こちらのソースを使わせてもらいました.
実際にモンテカルロ法で計算してみた
ソースコードはこんな感じに
using System; using jp.takel.PseudoRandom; namespace MCMethod { class Program { static void Main(string[] args) { //メルセンヌツイスター法で乱数発生 MersenneTwister mt = new MersenneTwister(); //実行時間を測る System.Diagnostics.Stopwatch sw = new System.Diagnostics.Stopwatch(); long time; //座標 double x; double y; //円の内部の点の数 double hit = 0; //点をプロットさせた回数 double trialNum = 0; //モンテカルロ法で計算した円周率 double pi; sw.Start(); for (long i = 0; i < 500000000; i++) { //乱数の最初の1000回分は使わない if (i > 1000) { x = mt.NextDouble(); y = mt.NextDouble(); trialNum++; if (Math.Sqrt(Math.Pow(x, 2) + Math.Pow(y, 2)) <= 1) { hit++; } } } pi = 4 * hit / trialNum; sw.Stop(); time = sw.ElapsedMilliseconds; // 表示 System.Console.WriteLine(pi); System.Console.WriteLine(pi - Math.PI); System.Console.WriteLine("実行時間: " + time + " ミリ秒"); System.Console.ReadLine(); } } }
