ScalaでProject Euler（100）

Scala

Problem 66

x - Dy2 = 1

これをPell方程式と呼びます。これには必ず整数解があることが知られています。この整数解が大きくなると、

x/y√(1-x-2) = √D

よりx/yは√Dの近づくことになります。√Dの連分数展開と同じですね。実際、Pell方程式の解は連分数展開を途中で打ち切ったものになります。例えば√3を連分数展開していくと、

1
1 + 1 / 1 = 2
1 + 1 / (1 + 1 / 2) = 5 / 3
…

ですが、2の分子をx、分母をyとすると、Pell方程式が成り立つことがわかります。このように、連分数展開していくといつか初期解が求まります。
この解は場合によって非常に大きくなります。ここでは本質がわかるようにBigIntを使いましょう。
ところで、Scalaでは

println (BigInt(2) * 3)

はOKでも

println (3 * BigInt(2))

はエラーなんですね。

import scala.math._

def gcd(m :Int, n :Int) :Int = if(n == 0) m else gcd(n, m % n)

def is_square(n :Int) = {
    val m = sqrt(n.toDouble).toInt
    m * m == n
}

def min_solution(D :Int) :BigInt = {
    val n = sqrt(D.toDouble).toInt
    val x0 = (1, n, 1, 0)
    
    def next(x :(Int,Int,Int,Int)) :(Int,Int,Int,Int) = {
        val (a, b, c, m) = x
        val m1 = (a * n + b) / c
        val e = m1 * c - b
        val (a1, b1, c1) = (a * c, c * e, a * a * D - e * e)
        val d = gcd(gcd(a1, b1), c1)
        (a1 / d, b1 / d, c1 / d, m1)
    }
    
    def equalto(x :(Int,Int,Int,Int)) =
        x._1 == x0._1 && x._2 == x0._2 && x._3 == x0._3
    
    def next2(x :(BigInt,BigInt,BigInt,BigInt,Iterator[Int])) = {
        val (a2, b2, a1, b1, it) = x
        val n = it.next
        (a2 * n + a1, b2 * n + b1, a2, b2, it)
    }
    
    def is_satisfied(x :(BigInt,BigInt)) =
        x._1 * x._1 - x._2 * x._2 * D == 1
    
    val it = Iterator.iterate(next(x0))(next).drop(1).map(_._4)
    val y0 = (BigInt(n), BigInt(1), BigInt(1), BigInt(0), it)
    val it2 = Iterator.iterate(y0)(next2).map(x => (x._1, x._2))
    it2.filter(is_satisfied).next._1
}

val M = 1000
println ((2 to M).filter(D => !is_square(D)).
                    map(D => (min_solution(D), D)).max._2)