モンテカルロでπを計算 part2
乱数でπを求めよう
では、点をプロットしていきましょう。まず、プロットする点の精度ですが、0.000～1.000とします。例えば、ある点のx軸が0.456でy軸が0.654なら以下の場所にプロットされます。
^y 1
|
-
|
|
- *
|
|
-
|
|
-
|
| x
|o 1
+-------|-------|-------|-------|------>
x軸とy軸に区切りを入れたので、大体の位置関係が分かると思います。精度は1000分の1ですが、テキスト文字での描画なので目安で見てもらうしかありません。
では、もっとたくさんの点をプロットしてみます。でたらめな点にするにはサイコロを振るなどがありますが、大量のでたらめ数字を作るのは大変です。そこでrpnの乱数を作る関数で、でたらめな数を作ります。
0.839
これで、精度が1000分の1のでたらめ数が得られます。これを2回繰り返せば、でたらめなx軸とy軸の点をプロットできます。まず、プロット数を10個にして実行してみましょう。
0.839 0.759
0.114 0.516
0.052 0.628
0.011 0.42
0.213 0.087
0.75 0.768
0.085 0.061
0.226 0.544
0.09 0.184
0.138 0.567
10個のでたらめな点が作られました。最初のプロットはx軸が0.839でy軸が0.759の位置ですね。同様にして残り9個をプロットしたグラフが次です。
^y 1
|
|
| * *
|
| *
| * *
| *
*
|
|
|
| *
| * x
|o * 1
+-------------------------------------->
rpnが選んだ、でたらめな数から作った10個の点がプロットしてあります。でたらめな割には何か塊があるような気がしますが、数が少ないとこのように規則があるように錯覚してしまいますね。
では、このグラフと上記の4分の1の円をnpdを使って重ねてみます。
| + +
| +
| + * *
| +
| * +
| * * +
| *
* +
| +
| +
|
| * +
| * +x
|o * 1
+-------------------------------------+>
見た目で確認すると、2つが円外で8つが円内ですね。rpnで計算してみます。
>rpn . * x . * + r <tmp | rpn 0 1 -c lookup | rpn -c count
8
rpnの計算でも、円の内側にある点の数は8つになりました。すると以下の式が成り立ちますね。
--- ＝ ---- ⇒ π ＝ 8÷10×4
4 10
早速、計算してみましょう。
3.2
3.2になりました。どうです。たった10個のでたらめな点から計算したπの値が本当のπの値(3.14)に相当近くなっていませんか。じゃあ、次は100個プロットしてみましょう。もっとπの精度が上がるでしょうか。グラフにしてみます。
^y*1 * ** * * * *
| * * **
| * * * **
| * * * * *
| * * *
| * * * ** *
|* ** ** * * * *
| ** ** * * * * **
* * ** *
| * * *
| * * * * * * *
| * * * * **
|* * * * * * *** *
| ** * * * * x
|o * * *** * 1
+-------------------------------------->
流石に100個となるとばらばらな感じにプロットされていますね。npdで4分の1の円と重ね合わせてみましょう。
| * + + * **
| * * * *+
| * * + * * *
| * * * +
| * * * **+ *
|* ** ** * * * + *
| ** ** * * * * **
* * ** + *
| * * * +
| * * * * * * +
| * * * * **
|* * * * * * *** *+
| ** * * * * +x
|o * * *** * 1
+-------------------------------------+>
目視で確認できる量ではありませんね。rpnで計算してみましょう。
>rpn . * x . * + r <tmp | rpn 0 1 -c lookup | rpn -c count
81
円内に81個がプロットされました。式に置き換えてみましょう。
--- ＝ ---- ⇒ π ＝ 81÷100×4
4 100
さあ、πはいくらになるでしょうか。rpnで計算です。
3.24
3.24になりました。10点のプロットよりも悪くなりましたね。でたらめな数なのでこんなこともあります。もっと数を増やせばより正確になるはずです。
どんどんプロット点を増やしてπへ近似
1000点をプロットしてみましょう。グラフは省略して、rpnの式だけを載せます。
>rpn 0 1 -c lookup <tmp | rpn -c count | rpn 1000 / 4 *
3.188
3.188です。やっと、3.1台にきました。次は1万点でやってみましょう。
>rpn 0 1 -c lookup <tmp | rpn -c count | rpn 10000 / 4 *
3.1604
3.1604です。また近づきました。では、5万だとどうでしょう。
>rpn 0 1 -c lookup <tmp | rpn -c count | rpn 50000 / 4 *
3.15976
3.15台まできました。でも、5万回プロットしてもアルキメデスの精度(3.14)には届きません。この方法では、道のりは遠いという感じですね。
番外編：コンピュータのπ計算についての補足
実際に、コンピュータでπを計算するときは、以下のような数式を使います。
π ＝ Σ { ------------------ } × 6
n=0 4n+1 2
2 ・n! ・(2n+1)
rpnで計算してみましょう。以下がそのrpn式ですが、長いので2行に分けて表示してあります。
===(この行の1行下からコピー)===
@s 2 @n * ! 2 4 @n * 1 + p @n ! . * * 2 @n * 1 + * / + #s \
@n 1 + #n @n @s 6 * -r 10
===(この行の1行上までコピー)===
上記のrpn式を一旦、テキストファイルに格納します。テキストファイルに付けるファイル名は何でもいいのですが、ここでは"pi.rpn"にしておきましょう。
では、さっそく計算です。数式からπはnが0から無限大までの総和になるのですが、とりあえずnは10までとします。以下のようにrpn式を格納したファイルを読み込ませることで計算できます。
1 3
2 3.125
3 3.13906
4 3.14116
5 3.14151
6 3.14158
7 3.14159
8 3.14159
9 3.14159
10 3.14159
πの計算式を10回繰り返しながら毎回、計算結果をレジスタに累積していきます。計算経過を見ると、n=7で単精度の限界に達しています。相当に速い収束ですね。
rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。
lookupはユーティリティパッケージに同梱されています。xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。
/