ビュフォンの針からπを近似する問題の面白かったとこ(C言語)

ビュフォンの針

ビュフォンの針 - Wikipedia

※適当なこと書いてる可能性が大いにあるのでお気をつけください．
平行線を一定間隔dで書いてそこに長さlの針を落とすと何本平行線に交わるかって問題．
授業でやって面白かったのでまとめてみたくなった．
l >= d　のときはなんか難しそうなので　
l < d　のときの解析解 P = 2l / πd　を使って針が平行線に交わる確率Pをπを使わずにモンテカルロ法で求める．そのあとに上の式を変形してπについて解くことでπの近似ができるって話．

さっそくできたコードを張ってみる

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

int main(){    
    int N, i, l = 3, d = 5, count = 0, flag_sin = 0;
    double x, theta, x0, y0, r, sin, P, pi;

    srand((unsigned int)time(NULL)); 
    printf("How many attempts:"); scanf("%d", &N);
    
    for(i = 0; i < N; i++){        
        x = ((double)d / 2) * ((double)rand() / RAND_MAX);         
        while(flag_sin == 0){            
            x0 = ((double)rand() / RAND_MAX); 
            y0 = ((double)rand() / RAND_MAX);
            r = sqrt(x0 * x0 + y0 * y0);            
            if(x0 * x0 + y0 * y0 <= 1){                 
                sin = y0 / r;
                flag_sin = 1;
            }
        }        
        flag_sin = 0;       
        if (x < (((double)l / 2) * sin)){
            count++;
        }
    }    
    P = (double)count / N;    
    pi = (2 * (double)l) / (P * d);

    printf("l = %d, d = %d, N = %d, count = %d\n", l, d, N, count);
    printf("P = %lf, pi = %lf\n", P, pi);

    return 0;
}

実行結果

1)
How many attempts:100000
l = 4, d = 20, N = 100000, count = 12746
P = 0.127460, pi = 3.138239

2)
How many attempts:100000
l = 3, d = 100, N = 100000, count = 1888
P = 0.018880, pi = 3.177966

3)
How many attempts:100000
l = 3, d = 100000, N = 100000, count = 2
P = 0.000020, pi = 3.000000

4)
How many attempts:100000
l = 3, d = 5, N = 100000, count = 37990
P = 0.379900, pi = 3.158726

5)
How many attempts:10000000
l = 3, d = 5, N = 10000000, count = 3821079
P = 0.382108, pi = 3.140474

なにがそんなに面白かったのか

針が平行線に交わっているかを判定するには

上の図のように針の中心から
一番近い変更線までの距離x
と
半分の針のy軸方向の長さ(ℓ / 2) * sinθ
の大小関係が
x < (ℓ / 2) * sinθ
になっていれば針が平行線と交わると言える．

これを使えば座標を使わなくても針が交わっているかどうか判定できるため，あとは
1. xを(0 <= x <= d / 2)の範囲でランダムに生成
2. θを(0 <= θ <= π/2)の範囲でランダムに生成
したら求めたい確率が計算できる．

1のxを(0 <= x <= d / 2)の範囲でランダムに生成は
x = (d / 2) * ((double)rand() / RAND_MAX);
で簡単にできる．
が，θの場合はそう簡単にいかない．

πを使わずにθを(0 <= θ <= π/2)の範囲でランダムに("均一に")生成するには
x = (π / 2) * ((double)rand() / RAND_MAX);
などとしないといけないが，これをするにはπの値が必要になる．
そのため，下の図のようにx, yを(0 <= x, y <= 1)の範囲でランダムに生成し，そのx, yからrを計算してx < (ℓ / 2) * sinθに必要なsinθをランダムに生成する．

ここで注意しなければならないのはsinθを計算していいのは点が円周内にある場合のみであるということ．
なぜかというとx, yを(0 <= x, y <= 1)の範囲でランダムに生成すると，
θ = 0, π / 2付近と
θ = π / 4付近
で点が少なくあらわれる範囲と点が多くあらわれる範囲があるということ．
つまりx, yを(0 <= x, y <= 1)の範囲でランダムに生成するとθは("均一に")ランダムに生成されないということ．

実際に，
一番現れやすい θ = π / 4 では 単位円の半径を1として√2の範囲があり，
一番現れにくい θ = 0, π / 2 では1の範囲しかない．
(ちなみに同じような理由で
sin = ((double)rand() / RAND_MAX); (0 <= sinθ <= 1)
などともできない)(たぶん)

(後述するが，これだとx1, y1はsinθをランダム生成するために使ってはいけない．)
ちょっと話が感覚的すぎるのでグラフをつくってみた．

X[i] to X[i + 1]でプロットされた点が"均一"にあるとしたら，Y[i] to Y[i + 1]でも同じだけ点あるはずだが，Xが一定なのに対し，Yは一定ではない．
つまり，x, yを(0 <= x, y <= 1)の範囲でランダムに生成すると角度θは
振り当てられやすい場所(上の画像だと0.38 / 0.38 + 0.32 + 0.22 + 0.08+…で振り当てられるなど)
と
振り当てられにくい場所(上の画像だと0.08 / 0.38 + 0.32 + 0.22 + 0.08+…で振り当てられるなど)
がある．

実は前のいきなりsinθをグリッド表示しだしたのはこれが言いたかったから

つまり範囲が半径の1で均一になる円周内で点をプロットしたときだけθが均一にランダム生成されたと言え，その時のsinθを使って求めたい確率Pを計算することができる．

あとは1つのランダムなsinθが生成されるまでx, yをランダムに生成するのをwhile文で書いて
P = count / N (ラプラスの着想)
で近似し，そこからπを近似した．