モンテカルロでπを求める part2 | 数学アラカルト [実践] | 逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓rpnの実践ウェブ   
rpn hacks! アールピーエヌ・ハックスサイトマップ

rpn | 実戦 | 数学アラカルト | 数学に関する話題をrpnで探求!数字と遊んでみよう。

HOME > 実践 > 数学アラカルト > モンテカルロでπを求める part2 hatena yahoo buzzurl livedoor del.icio.us nifty newsing twitter facebook rss ソーシャルブックマーク

モンテカルロでπを求める part2

前のページに戻るLinkIcon

乱数でπを求めよう

 では、点をプロットしていきましょう。まず、プロットする点の精度ですが、0.000~1.000とします。例えば、ある点のx軸が0.456でy軸が0.654なら以下の場所にプロットされます。

  >rpn .456 .654 | xyp -x,1 -y,1 -s.2,.2 -m
  ^y 1
  |
  -
  |
  |
  -                *
  |
  |
  -
  |
  |
  -
  |
  |                                      x
  |o                                     1
  +-------|-------|-------|-------|------>


x軸とy軸に区切りを入れたので、大体の位置関係が分かると思います。精度は1000分の1ですが、テキスト文字での描画なので目安で見てもらうしかありません。

では、もっとたくさんの点をプロットしてみます。でたらめな点にするにはサイコロを振るなどがありますが、大量のでたらめ数字を作るのは大変です。そこでrpnの乱数を作る関数で、でたらめな数を作ります。

  >rpn 1000 ? 1000 /
  0.839


これで、精度が1000分の1のでたらめ数が得られます。これを2回繰り返せば、でたらめなx軸とy軸の点をプロットできます。まず、プロット数を10個にして実行してみましょう。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 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個をプロットしたグラフが次です。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 10 | xyp -x,1 -y,1 -m
  ^y 1
  |
  |
  |                            *  *
  |
  | *
  |    *  *
  |   *
  *
  |
  |
  |
  |  *
  |       *                              x
  |o *                                   1
  +-------------------------------------->


rpnが選んだ、でたらめな数から作った10個の点がプロットしてあります。でたらめな割には何か塊があるような気がしますが、数が少ないとこのように規則があるように錯覚してしまいますね。

では、このグラフと上記の4分の1の円をnpdを使って重ねてみます。

  ^y +   +   +
  |              +  +
  |                     +
  |                        +   *  *
  |                          +
  | *                           +
  |    *  *                       +
  |   *
  *                                 +
  |                                   +
  |                                    +
  |
  |  *                                  +
  |       *                             +x
  |o *                                   1
  +-------------------------------------+>


見た目で確認すると、2つが円外で8つが円内ですね。rpnで計算してみます。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 10 >tmp
  >rpn . * x . * + r <tmp | rpn 0 1 -c lookup | rpn -c count
  8


rpnの計算でも、円の内側にある点の数は8つになりました。すると以下の式が成り立ちますね。

   π      8
  --- = ----    ⇒   π = 8÷10×4
   4      10


早速、計算してみましょう。

  >rpn 8 10 / 4 *
  3.2


3.2になりました。どうです。たった10個のでたらめな点から計算したπの値が本当のπの値(3.14)に相当近くなっていませんか。じゃあ、次は100個プロットしてみましょう。もっとπの精度が上がるでしょうか。グラフにしてみます。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 100 | xyp -x,1 -y,1 -m
  ^y*1  *          ** *  *     *       *
  |            *               *      **
  |   *          * *   **
  |   *    *                   *  * *
  |             *   * *
  | *               *   *     **    *
  |*  **  **      *   * *           *
  |  **  **    * *    *    * **
  *     *          **                  *
  |   *     *                  *
  |               *  * *  *     *    * *
  | *             *       * * **
  |* *   *   *  *     *  ***           *
  |      **      *     *  *        *     x
  |o *  *     ***        *               1
  +-------------------------------------->


流石に100個となるとばらばらな感じにプロットされていますね。npdで4分の1の円と重ね合わせてみましょう。

  ^y*+  *+   +     ** *  *     *       *
  |            * +  +          *      **
  |   *          * *   *+
  |   *    *               +   *  * *
  |             *   * *      +
  | *               *   *     **+   *
  |*  **  **      *   * *         + *
  |  **  **    * *    *    * **
  *     *          **               +  *
  |   *     *                  *      +
  |               *  * *  *     *    * +
  | *             *       * * **
  |* *   *   *  *     *  ***           *+
  |      **      *     *  *        *    +x
  |o *  *     ***        *               1
  +-------------------------------------+>


目視で確認できる量ではありませんね。rpnで計算してみましょう。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 100 >tmp
  >rpn . * x . * + r <tmp | rpn 0 1 -c lookup | rpn -c count
  81


円内に81個がプロットされました。式に置き換えてみましょう。

   π     81
  --- = ----    ⇒   π = 81÷100×4
   4     100


さあ、πはいくらになるでしょうか。rpnで計算です。

  >rpn 81 100 / 4 *
  3.24


3.24になりました。10点のプロットよりも悪くなりましたね。でたらめな数なのでこんなこともあります。もっと数を増やせばより正確になるはずです。

どんどんプロット点を増やしてπへ近似

 1000点をプロットしてみましょう。グラフは省略して、rpnの式だけを載せます。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 1000 | rpn . * x . * + r >tmp
  >rpn 0 1 -c lookup <tmp | rpn -c count | rpn 1000 / 4 *
  3.188


3.188です。やっと、3.1台にきました。次は1万点でやってみましょう。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 10000 | rpn . * x . * + r >tmp
  >rpn 0 1 -c lookup <tmp | rpn -c count | rpn 10000 / 4 *
  3.1604


3.1604です。また近づきました。では、5万だとどうでしょう。

  >rpn 1000 ? 1000 / 1000 ? 1000 / -r 50000 | rpn . * x . * + r >tmp
  >rpn 0 1 -c lookup <tmp | rpn -c count | rpn 50000 / 4 *
  3.15976


3.15台まできました。でも、5万回プロットしてもアルキメデスの精度(3.14)には届きません。この方法では、道のりは遠いという感じですね。

番外編:コンピュータのπ計算についての補足

 実際に、コンピュータでπを計算するときは、以下のような数式を使います。

        ∞          (2n)!
  π = Σ { ------------------ } × 6
        n=0    4n+1   2
              2   ・n! ・(2n+1)


rpnで計算してみましょう。以下がそのrpn式ですが、長いので2行に分けて表示してあります。

 【π計算のrpnプログラム(pi.rpn)】
  ===(この行の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式を格納したファイルを読み込ませることで計算できます。

  >rpn -r 10 <pi.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で単精度の限界に達しています。相当に速い収束ですね。

実践数学アラカルトに戻るLinkIcon

警告rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。

情報lookupはユーティリティパッケージに同梱されています。xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。

数学アラカルト

数字と遊んでみよう

※実践コーナーのTOP

紹介 rpnの利用シーンはこちら…

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

part2

数字と遊んでみよう

※実践コーナーのTOP

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

書籍紹介

記事に関連した書籍

本ウェブサイトで扱った話題に関連した書物で、スタッフが実際に読了したものを紹介。

書籍数学の書籍
数の世界は思ったよりもエキサイティング。

  • 書籍統計の書籍
  • ビジネスで統計が使えるととっても有利。

書籍投資の書籍
失敗しない投資には広範囲で実践的な知識が必要。

警告バックスラッシュはエンマークに読み替えてください( IEのみ)。
バックスラッシュとエンマーク

警告文字で作られた図表や式が崩れることがあります。ブラウザによっては固定幅フォントをMSゴシックにするときれいに表示されます。それでも崩れる場合は図表や式をメモ帳にコピー後、閲覧下さい。

警告rpn試用版と標準版(2kリビジョン)はダブルクォートで囲ってください。

rpn 1 2 + ⇒ rpn "1 2 +"
rpn 1 -c foo ⇒ rpn "1" -c "foo"

ダブルクォートは省略できることが多いのですが、慣れない間は囲んだほうが無難です。なお、本ウェブサイトの記事ではrpn標準版(98リビジョン)を使用しているため囲っていません。詳しくは技術サポートの「rpn TIPS参照ください。

注意rpnの障害情報はこちら

警告rpn試用版の場合、複雑なプログラムや処理時間のかかるプログラムの一部には動作しないものがあるかもしれません。あくまで無料提供であることを勘案・了承ください。rpn標準版は、すべてのプログラムが動作します。