目次
Kotlin で、 モンテカルロ法を使って円周率を計算しました。
モンテカルロ法
乱数を用いて近似値を計算する方法です。
円周率の計算イメージ
平面に、原点を中心とする単位円(原点を中心とする半径1の円)があるとします。 4点 ( (x, y) = (1, 1), (1, -1), (-1, -1), (1, -1) ) を順に直線で結んでできる正方形の中にランダムに点を描画し、その点のうち単位円の中に描画された点の割合から円周率を計算します。
環境
- Kotlin 1.1.2-2
コード
次のようにしました。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
package com.improve_future.math object Main { @JvmStatic fun main(args: Array<String>) { val unitCircle = Circle(0.0, 0.0, 1.0) val m = 1000000L var c = 0L for (l: Long in 1L..m) { val p = Point2d(random(), random()) if (unitCircle.contain(p)) c += 1L } println(c * 4.0 / m) } fun random(): Double { return Math.random() * 2 - 1 } } |
|
1 2 3 4 5 6 7 |
package com.improve_future.math class Circle(val x: Double, val y: Double, val r: Double) { fun contain(point: Point2d): Boolean { return Math.pow(point.x - x, 2.0) + Math.pow(point.y - y, 2.0) < Math.pow(r, 2.0) } } |
|
1 2 3 |
package com.improve_future.math data class Point2d(val x: Double, val y: Double) |
円を表現するクラス Circle を作っていますが、 単位円にしか使いません。 点を表現するクラス Point2d も作りましたが、 クラスにしなくても計算はできます。
まず単位円のオブジェクトを作ります。 そして ( -1 leq x lt 1, -1 leq y lt 1 ) となるように点をランダムに作成して、 そのうち単位円の中に存在する点の数を計算します。 理論上、 単位円の面積は ( pi ) になり、 点が存在しうる範囲の面積は ( 4 ) になります。 よって、 ( 4 ) に 単位円の中の点の割合を掛けて円周率を計算します。
計算例
m = 1000000000L として実行したところ、 次のように表示されました。
|
1 |
3.141752604 |
生成される乱数に偏りがあると、実際の値との乖離が大きくなります。
格子点でやってみる
なにも乱数にしなくても、格子点でやってみても円周率は計算できそうですね。 main関数を変更して実行してみます。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
package com.improve_future.math object Main { @JvmStatic fun main(args: Array<String>) { val unitCircle = Circle(0.0, 0.0, 1.0) val m = 35000L println(m) var c = 0L for (a: Long in 1L..m) { for (b: Long in 1L..m) { val p = Point2d(2.0 * a / m - 1, 2.0 * b / m - 1) if (unitCircle.contain(p)) c += 1L } } println(c * 4.0 / m / m) } } |
これは乱数を使っていませんから、何度やっても同じ値が出ます。 しかもモンテカルロ法よりも円周率に近い値が出ました。
実行結果
|
1 |
3.1415914351020406 |
