ScalaでProject Euler（56）

Problem 31

さらに速い方法があります。

P_k(x) = 1 / (1 - x^k)

なので、

(1 - x)(1 - x²)...(1 - x¹⁰⁰)(1 - x²⁰⁰)

を計算して逆数を取るだけです。本当は証明が必要なのですが手抜きします。
N = 5000で前回の方法より2桁速かったです。O(N²)がO(N)ですからね。

import scala.testing.Benchmark

type Poly = Array[Long]

val coins = List(1, 2, 5, 10, 20, 50, 100, 200)
val N = 5000
val M = coins.sum

def unit() : Poly =
    Iterator.range(0, N + 1).map(k => if(k == 0) 1L else 0L).toArray

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(N + 1)
    val a = f.clone
    for(k <- 0 to N) {
        val c = a(k) / g(0)
        h(k) = c
        for(l <- 0 to M.min(N - k))
            a(k + l) -= c * g(l)
    }
    h
}

def solve() = {
    val f = coins.foldLeft(unit)(mul)
    val g = div(unit, f)
    println (g(N))
}

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

println (Test.runBenchmark(10))