ScalaでProject Euler（101）

Problem 66

Pell方程式の最小解と連分数には驚くべき関係があります。たとえば、Problem 64で例に挙げられた√23は、

√23 + 4 = [ 8; 1, 3, 1, 8, 1, 3, 1 ... ]

ですが、1周期で打ち切ったとき、

[ 4; 1, 3, 1 ] = 24 / 5

で、この(24, 5)が

24² - 23 * 5² = 1

で初期解になっています。別の例でD = 61だと、

√61 + 7 = [ 14, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14, ... ]

で、1周期で打ち切ったとき、29718 / 3805となり、

29718² - 61 * 3805² = -1

となります。実は、連分数を1周期で打ち切ったとき、必ず1か-1になります。1ならばそれが初期解で、-1なら

a² - Db² = -1
(a + b√D)(a - b√D) = -1
(a + b√D)²(a - b√D)² = 1
(a² + Db² + 2ab√D)(a² + Db² - 2ab√D) = 1

(a² + Db², 2ab)が初期解となります。

また、連分数の周期は、最初の14がもう一度出てくる手前までなので、この数だけをチェックすればよいです。
前回より4.5倍ほど速くなりました。Pythonでは2倍程度でした。Scalaは多倍長整数が遅いからではないかと思われます。

import scala.math._
import scala.testing.Benchmark

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 m = sqrt(D.toDouble).toInt
    
    def next_n(x :(Int,Int,Int,Int)) :(Int,Int,Int,Int) = {
        val (a, b, c, n) = x
        val n1 = (a * m + b) / c
        val e = n1 * 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, n1)
    }
    
    def next_frac(x :(BigInt,BigInt,BigInt,BigInt), n :Int) = {
        val (a2, b2, a1, b1) = x
        (a2 * n + a1, b2 * n + b1, a2, b2)
    }
    
    val x0 = (1, -m, 1, 0)
    val it = Iterator.iterate(x0)(next_n).map(_._4).takeWhile(m * 2 !=)
    val f0 = (BigInt(m), BigInt(1), BigInt(1), BigInt(0))
    val f = it.foldLeft(f0)(next_frac)
    val (x, y) = (f._1, f._2)
    if(x * x - y * y * D == 1)
        x
    else
        x * x + y * y * D
}

def solve() = {
    val M = 100000
    println ((2 to M).filter(D => !is_square(D)).
                    map(D => (min_solution(D), D)).max._2)
}

object Test extends Benchmark {
    def run() = solve
}

println (Test.runBenchmark(10))