対象を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)