pi で pi を求めてみた。(Scala自習会9日目)
閉じる
閉じる

新しい記事を投稿しました。シェアして読者に伝えましょう

×

pi で pi を求めてみた。(Scala自習会9日目)

2015-12-23 17:15
    ○学生時代に、ゆとり世代ならおなじみの「総合的な学習の時間」というのがあって、その時に数学の先生から、乱数で円周率を求める事が出来る、モンテカルロ法の存在を教わりました。

    それ以降、取りあえず新しい言語を覚えたら円周率を求めてみるということをやっているので、Scalaでもやってみました。

    で、完成したソースがこちら。元となる円周率は、100万桁を公開していらっしゃるサイトがあり、そこから直接拝借することにしています。(下記ソース参照 gistはこちら

    import skinny.http._

    object Main extends App {
    val response: Response = HTTP.get("http://www.geocities.jp/f9305710/PAI1000000.html")
    var r = response.textBody
    r = r.substring(r.indexOfSlice("141592"), r.lastIndexOfSlice("</FONT>")).withFilter(_.isDigit).map(c => c)

    var incount = 0;
    var oncount = 0;
    var outcount = 0

    (r grouped 32) foreach { s =>
    if (s.length == 32) {
    print("" + s.substring(0, 15) + " " + s.substring(16, 31) + " ")
    var p: Double = Math.pow(("0." + s.substring(0, 15)).toDouble, 2) + Math.pow(("0." + s.substring(16, 31)).toDouble, 2)
    print(p)
    p match {
    case q if q < 1 => incount += 1
    case 1 => oncount += 1
    case q if q > 1 => outcount += 1
    }
    println(" " + incount + " " + oncount + " " + outcount + " " + (4.0 * incount / (incount + outcount)))
    }
    }
    }


    実際にでた数字は、

    3.1425646275915025

    下2桁しか合ってない。でも0.03%の差異なら、さほど悪くない?プロットしてみても値が偏っているようにも見えないし。なかなか難しいですね。








    そろそろ、「え、円周率を100万桁使ってて、なんで円周率2桁しか出せないの?」って突っ込まれそうですが。

    (ここから私の数十年前のうろ覚えを元に説明。ちゃんと知りたい人はググってください。)

    モンテカルロ法というのは、乱数を使って複雑な形の面積を求めたり、簡単に計算では求められない数字を概算するための方法らしいです。



    たまたま適当に12個の点を上の正方形の中に置いたとします。
    すると、その中の円の中に偶然入った9点と、その中の円に偶然入らなかった3点が出てきます。
    この点の数の比率を計算することで、正方形の中の円の面積を大体計算出来るはずと考えられるのです。
    上の円の面積を、pi * r^2 とすると、正方形の面積は、(2r)^2となります。
    pi * r^2 : (2r)^2 = 9 : 9+3
    計算すると、
    pi = 3
    となり、円周率は3に近い数字であることがわかります。

    もちろん、点の数を増やせば増やすほど、精度は上がると考えられるので、その点の数を増やすには、多くの乱数が必要になって、多くの乱数といえば、、、円周率って、かなり理想的な乱数と聞いたことがある。。。。よし、円周率を使って円周率を求めよう!

    (この時点で本末転倒なのは置いておいて。)

    上のプログラムでは、
    1. 円周率100万桁をとあるホームページから拝借してくる(数字以外の文字列を排除して数字の文字列に整形)
    2. 32文字ずつに分割して、更にそれを半分にして、それぞれ、0.[1-16桁目] と、 0.[17-32桁目] という0以上1未満の数字に変換して、(x, y) 座標を作成
    3. r = 1として、 円の公式、x^2 + y^2 = 1を利用して、 x^2 + y^2 が、1より大きいか、1より小さいかで、作成した(x, y) 座標が円の外か、円の中かを判定
    4. それを数字が無くなるまで繰り返す
    ということをしています。
    16桁というのは、Doubleの大体の有効桁数ということで。
    100万/32 = 大体3万サンプルほど使っても、下2桁しか正格に出せないということで、あまり使い物になりませんし、もっと早く正確な数字を出す方法は一杯あるので、実用性はないですが、○学生でも分かる仕組みで、簡単にプログラムが書けるという意味で私は気に入ってます。
    ちなみに、昔試したときは、普通のC言語とかの乱数を使うと下一桁ぐらいが限界で、Excelだと突拍子もない数字になりました。今はメルセンヌツイスターとかいうカッコイイ物があるのでもっと精度出るのかも知れないです。

    そして、Raspberry Pi には、物理乱数の生成モジュールが内蔵されているらしいので、「piでpiを求めてみた」の第二弾はもう決まりですね。
    広告
    コメントを書く
    コメントをするには、
    ログインして下さい。