ガウス=ルジャンドルのアルゴリズムとは
円周率を計算するアルゴリズムです。
ソースコード
今回はすばやく実装することを目指したので、Double
で実装しています。故に精度は有りません。
計算はアルゴリズムそのままな書き方をしているので解説しません。
/**
* 計算用の入れ物
*/
private data class CalcData(
val a: Double,
val b: Double,
val t: Double,
val p: Int
)
/**
* 2乗
*/
private fun sqr(d: Double) = d * d
/**
* 二乗根
*/
private fun sqrt(d: Double) = Math.pow(d, 0.5)
/**
* 計算の本体
*/
fun gaussLegendre(n: Int) = generateSequence(
// 初期値
CalcData(1.0, 1.0 / sqrt(2.0), 0.25, 1),
// 更新関数
{(a, b, t, p) ->
val aNext = (a + b) / 2.0
val bNext = sqrt(a * b)
val tNext = t - sqr(a - aNext) * p.toDouble()
val pNext = p.shl(1)
CalcData(aNext, bNext, tNext, pNext)
}
).take(n).last().let { sqr(it.a + it.b) / (it.t * 4.0) }
実行結果
大体n=4
でDouble
の精度の限界に達します。
テストコード
import org.junit.jupiter.api.Test
class TestGaussLegendle {
@Test
fun testGaussLegendle() {
println(gaussLegendre(3)) // 3.141592646213543
println(gaussLegendre(4)) // 3.141592653589794
println(gaussLegendre(5)) // 3.141592653589794
println(Math.PI) // 3.141592653589793
}
}