ScalaでProject Euler(110)

Problem 73

対象を1/3より大きく1/2より小さいより、1/3以下とするほうが考えやすそうなのでそうしましょう。1/2より小さいは、前問の1より小さいから1/2の1個を引いて半分にするだけです。

例えば、d = 40で考えましょう。取りうるnは1〜13です。ここから2の倍数を削ります。5の倍数も同じです。しかし、10の倍数は重複して削っているのでその分は足します。まとめると、まとめると分母40の分数の個数は、

[ 40/3 ] - [ 20/3 ] - [ 8/3 ] + [ 4/3 ]

これはメビウス関数μ(m)を用いて、

Σd|mμ(m)[ d/(3m) ]

と書けます。メビウス関数μ(m)は、mが平方数を因子に含んでいれば0、そうでなければ素数の個数が奇数なら-1、偶数なら1です。4や20は平方数を含んでいるので、40の約数でも上の式に現れないわけです。

これを100まで足し合わせてみます。まず1の倍数同士、2の倍数同士を足し合わせましょう。

1の倍数 : [ 1/3 ] + [ 2/3 ] + ... + [ 100/3 ]
2の倍数 : [ 1/3 ] + [ 2/3 ] + ... + [ 50/3 ]
3の倍数 : [ 1/3 ] + [ 2/3 ] + ... + [ 33/3 ]

33の倍数 : [ 1/3 ] + [ 2/3 ] + ... + [ 3/3 ]

そしてこれを符号に気を付けて足し合わせると、

f(n) = [ 1/3 ] + [ 2/3 ] + ... + [ n/3 ]

として、

Σ1≤n≤33μ(n)f(100/n)

となります。メビウス関数はエラトステネスのふるいと同等の計算量なので、O(NloglogN)となります。前回の解法はO(N2logN)程度でしょうか。


def div_pow(n :Int, d :Int) :Int =
    if(n % d != 0) n else div_pow(n / d, d)

def sieve(max_n :Int) :Array[Int] = {
    val a = new Array[Int](max_n + 1)
    for(n <- 1 to max_n) a(n) = n
    for(p <- Iterator.from(2).takeWhile(p => p * p <= max_n) if a(p) == p;
        n <- Iterator.range(p, max_n + 1, p))
            a(n) = p
    a
}

def phi(n :Int) :Int =
    if(n == 1)
        1
    else {
        val p = a(n)
        val m = div_pow(n, p)
        n / m / p * (p - 1) * phi(m)
    }

val N = 1000000
val a = sieve(N)
println (Iterator.range(2, N + 1).map(n => phi(n).toLong).sum)