ScalaでProject Euler（57）

Problem 31

もっと速くならないでしょうか。例えば一般項を求めるというのはどうでしょう。

Q₁(x) = P₁(x)

の一般項はa_k⁽¹⁾ = 1ですね。

Q₂(x) = P₁(x)P₂(x)

は、kが奇数なら、

a_k⁽²⁾ = (k + 1) / 2

は、kが偶数なら、

a_k⁽²⁾ = k / 2 + 1

となります。ここまでは簡単です。

Q₅(x) = P₁(x)P₂(x)P₅(x)

は、k = 10qなら、

a_k⁽⁵⁾ = a₀⁽²⁾ + a₅⁽²⁾ + ... + a_10q⁽²⁾
= (2 + 12 + ... + 10q + 2) / 2 + (6 + ... + 10q - 4) / 2
= 5q(q + 1) / 2 + q + 1 + 5q(q + 1) / 2 - 2q
= 5q² + 4q + 1

このような計算を計算機代数的に行っていって、Q₂₀₀(x)の一般項を出せばあとは式にNを代入するだけです。しかし、難しそうですよね。やれなくはないけど大変そう。

そうずっと思っていました。しかし、ひどい抜け道を思いつきました。3つの種類のコインなら一般項は2次式でした。ならば8つの種類のコインでは7次式でしょう。そして、200の剰余で一般項は分類されます。7次式ならばそれが8点あれば求まります。Nの200の剰余が0ならば、0, 200, 400, ... , 1400の8点を使って7次式を求めればよいでしょう。それはすでに求まっています。8点を通る7次方程式を求める方法は簡単です。まず、3点を通る2次式は、

です。3点を通ることが簡単に確認できると思います。一般にn点を通るn - 1次式は、

です。このxにNを代入すれば求める答えが出ます。N = 10⁸でも数ミリ秒で出ます。

import scala.testing.Benchmark

class Fraction(num :BigInt, den :BigInt) {
    private val (numerator, denominator) = init(num, den)
    
    def init(num :BigInt, den :BigInt) = {
        val d = gcd(num, den).abs
        if(den > 0)
            (num / d, den / d)
        else
            (-num / d, -den / d)
    }
    
    def *(f :Fraction) =
        new Fraction(numerator * f.numerator, denominator * f.denominator)
    
    def *(n :BigInt) =
        new Fraction(numerator * n, denominator)
    
    def +(f :Fraction) =
        new Fraction(numerator * f.denominator + denominator * f.numerator,
                                                denominator * f.denominator)
    
    override def toString() =
        if(denominator != 1)
            numerator.toString + "/" + denominator.toString
        else
            numerator.toString
    
    def gcd(m :BigInt, n :BigInt) :BigInt =
        if(n == 0) m else gcd(n, m % n)
}

type Poly = Array[Long]

val coins = List(1, 2, 5, 10, 20, 50, 100, 200)
val N = 1e8.toInt
val L = 200 * coins.size
val M = coins.sum

def unit() = {
    val a = new Poly(L + 1)
    a(0) = 1
    for(k <- 1 to L) a(k) = 0
    a
}

def mul(f :Poly, n :Int) = {
    val a = f.clone
    for(k <- Iterator.range(0, M - n + 1)) {
        a(k + n) -= f(k)
    }
    a
}

def div(f :Poly, g :Poly) = {
    val h = new Poly(L + 1)
    val a = f.clone
    for(k <- 0 to L) {
        val c = a(k) / g(0)
        h(k) = c
        for(l <- 0 to M.min(L - k))
            a(k + l) -= c * g(l)
    }
    h
}

def solve() = {
    val f = coins.foldLeft(unit)(mul)
    val g = div(unit, f)
    val xs = Array.range(N % 200, L, 200)
    val ys = xs.map(g(_))
    
    println ((0 to 7).map(k => (0 to 7).filter(k !=).
                map(l => new Fraction(N - xs(l), xs(k) - xs(l))).
                reduceLeft(_ * _) * ys(k)).
                reduceLeft(_ + _))
}

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

println (Test.runBenchmark(20))