TSperl 発言

極悪です。長文かも。ネタです。悩んでなんかいません（笑）

Perl の乱数は怪しいじゃないか、という話の中で、モンテカルロ
でπを求める話が出てきました。これは正方形内に打ったランダム
な点のうち、内接している円内にあるものの割合を求めれば円の面
積が分かってπも求まるという話。このとき乱数が一様でないと正
確な結果が出ないわけです。

僕は以前試したことがありますが、乱数の種の選び方に問題がある
のかぜんぜん精度が出ませんでした。いくら時間をかけても有効桁
が伸びず、小数点以下４〜５桁くらいでフラフラしてたと思います。
そのときは「なんだこりゃ、これなら乱数ではなく等間隔の点のほ
うがマシかな」と思いつつ諦めました。

んで、今回話が出たので再度やってみました。

計算は半径１の円（の第一象限）で考えます：

    $pi = 4 * $p / $a;
    
        $p : 円の内部に落ちた点の数
        $a : 総点数（正方形内に落とした点の数）



まずは普通にモンテカルロ：

srand 1;
my $a;      # whole area
my $p;      # inner area
my $acc = 0;
my $n   = 1;
while($acc < 8){
    $n  += $n;
    my $da  = $n*$n/4;
    my $dp;
    foreach(1..$da){
        my $x   = rand;
        my $y   = rand;
        ++$dp if $x*$x+$y*$y <= 1;
    }
    my $pi  = ($p += $dp)/($a += $da)*4;
    printf "(%d) %12.10lf\n",
                $acc = length $p,$pi;
}

無用に複雑なことをしてますが、それは下の等間隔メッシュ版を先
に作ったから。



次に乱数を使わず、等間隔にメッシュを切った場合（今回はこれを
試したかった）：

my $a;      # whole area
my $p;      # inner area
my $acc = 0;
my $n   = 1;
while($acc < 8){
    $n  += $n;
    my $rr  = $n*$n;
    my $da  = $rr/4;
    my $dp;
    for(my $ix = 1; $ix <= $n; $ix += 2){
        for(my $iy = 1; $iy <= $n; $iy += 2){
            ++$dp if $ix*$ix+$iy*$iy <= $rr;
        }
    }
    my $pi  = ($p += $dp)/($a += $da)*4;
    printf "(%d) %12.10lf \n",
                $acc = length $p,$pi;
}

メッシュを倍々で細かくしていったときのπの値を表示するように
なっています。やたら複雑になってる原因は、メッシュを倍にした
ときに追加分（$dp,$da）だけを計算するようにしているため。終
了条件が length $p < 8 になってるのは、こうしとけば有効桁８
桁のπが求まるんじゃないか、という朝知恵。



さて結果ですが、両方とも五十歩百歩でした。

モンテカルロでは
%timer "perl mont_pi.pl"
(1) 4.0000000000
〜中略〜
(7) 3.1414205516
(8) 3.1412784329
36.521000 sec

等間隔メッシュでは
%timer "perl mesh_pi.pl"
(1) 4.0000000000
〜中略〜
(7) 3.1416251148
(8) 3.1416038743
23.577000 sec

やはりモンテカルロよりも等間隔メッシュのほうが精度は良い (*1)
ようですが、贔屓目に見ても５桁、とても目標の８桁まではいきませ
ん。乱数生成ルーチンをモジュールの Math::Random::MT (*2) で置
き換えてみましたが、

%timer "perl montMT_pi.pl"
(1) 0.0000000000
〜中略〜
(7) 3.1418346847
(8) 3.1416358820
207.684000 sec
%

となって、これもどんぐりの背比べです（というか時間かかりす
ぎ）。経験上、一晩ほっといてもこれ以上精度は出ないと思います。
スピードはともかく、理屈は合ってるはずなのに なぜ精度が出な
いのでしょうか。




(*1) だいたいこういう問題に乱数を使うのが間違ってる気がする。
乱数ってのは規則性が邪魔になる場合に使うものじゃないかな。
逆に言えば、規則的に打った点の位置を変調してから使えばモンテ
カルロで乱数を使わなくても済むと思う。

(*2) ActivePerl なら ppm でインストールできます。
>The Mersenne Twister is a pseudorandom number generator
>developed byMakoto Matsumoto and Takuji Nishimura. It is
>described in their paper at
><URL:http://www.math.keio.ac.jp~nisimura/random/doc/mt.ps>.
だそうな。
-- 
FZH01112@..., http://www1.u-netsurf.ne.jp/~dune/Perl_82_C5pi.html?