モンテカルロ法で円周率を求める

d言語

「モンテカルロ法を実装しなきゃ！」という気分になった。

まずモンテカルロ法がなんなのかわからんかったので、調べるところから。適当にぐぐった。「確率で積分がどうしたこうしてほげほげ」。よくわからんな。円周率を求めることができるらしい。なんか聞いたことあるねそれ！ よくわからんまま適当にやってみよう。

マニュアルのstd.randomのページに「Cのrandよりランダム性高いよ！」と書いてあったので、Cのrandとphobosのrandを切り替えられるようにしてみた。

% repeat 5 (./pi 50000000 c)
testNum=50000000, mode=c,               pi=3.1419588000
testNum=50000000, mode=c,               pi=3.1416315200
testNum=50000000, mode=c,               pi=3.1417908000
testNum=50000000, mode=c,               pi=3.1414636800
testNum=50000000, mode=c,               pi=3.1417736000
% repeat 5 (./pi 50000000 phobos)
testNum=50000000, mode=phobos,        pi=3.1415154400
testNum=50000000, mode=phobos,        pi=3.1416636000
testNum=50000000, mode=phobos,        pi=3.1419204000
testNum=50000000, mode=phobos,        pi=3.1416398400
testNum=50000000, mode=phobos,        pi=3.1413786400

わーすげえや円周率っぽい！

うーん。っつーかこんなん、doubleとかintとかuintの長さとか、randの返す範囲とかそういうのも考慮に入れて考えなきゃいけないから、単純に「こっちのが円周率に近かったからよい擬似乱数！」って喜んでいいもん違うだろうし。以下ソース。

import std.stdio;
import std.string;
import std.random;
import std.c.stdlib;
import std.c.time;
double getPi(int testNum, double function() random)
{
  double pi = 0;
  int count = 0;
  for(int i=0;i<testNum;++i){
    double x = random();
    double y = random();
    if(x*x+y*y <= 1) ++count;
  }
  return count/cast(double)testNum*4.0;
}
double phobosRand(){
  double rand = std.random.rand();
  rand /= uint.max;
  return rand;
}
// stdlib.hに書いてる値もってくる。
const int RAND_MAX = 2147483647;
double cRand(){
  double rand = std.c.stdlib.rand();
  rand /= RAND_MAX;
  return rand;
}
double function() getRand(char[] type)
{
  double function() random;
  if(type=="phobos"){
    random = &phobosRand;
  }
  else if(type=="c"){
    srand(cast(uint)time(null));
    random = &cRand;
  }
  else{
    random = &phobosRand;
  }
  return random;
}
int main(char[][] args)
{
  if(args.length==1){
    fwritefln(stderr, "usage: %s testNum [phobos|c]", args[0]);
    return 1;
  }

  int n = cast(int)std.string.atoi(args[1]);
  char[] type = args.length==3 ? args[2] : "phobos";
  double function() random = getRand(type);

  writefln("testNum=%d, mode=%s,\t\tpi=%.10f", n, type, getPi(n, random));
  return 0;
}

べつにこれC言語でもほとんど同じに書けるなそういえば。あんましかっとんだ機能使わずにD言語で書くのは、C言語で書くプログラムのプロトタイプには良いのかもしれないなあ。

test