φ関数といえば、nの約数に亘って和を取るとnになる、という定理があります。Pythonで書くとこんなイメージです。
n == sum(map(phi, filter(lambda d: n % d == 0, range(1, n + 1))))
これを利用できないでしょうか。
1 + … + 100
にこの定理を適用すると、
1 + … + 100
= (φ(1)) + (φ(1) + φ(2)) + … + (φ(1) + … + φ(50) + φ(100))
= 100φ(1) + 50φ(2) + 33φ(3) + … + φ(100)
となります。
S(n) = φ(1) + … + φ(n)
とすると、
S(100) = 100 * 101 / 2 - φ(50) - … -49φ(2) - 99φ(1)
となり、φを100の半分の50まで求めればいいことになります。
# sieveとgen_phiは下のコードと同じ
N = 10 ** 6
M = N / 2
primes = sieve(M)
L = len(primes)
print N * (N + 1) / 2 - sum(imap(
lambda x: x[1] * (N / x[0] - 1), gen_phi(M))) - 1
これだけでもかなり速くなります。
これを再帰的に適用するとさらに速くなります。ただし、適当な大きさより小さいところでは個々のφを用いて計算します。この適当な大きさは√Nより少し大きいくらいがいいようです。さらにメモ化を使うと速くなります。
from itertools import ifilter, imap, takewhile
from math import sqrtdef sieve(max_n):
L = (max_n + 1) / 2
a = [ True ] * L
for k in ifilter(lambda p: a[p], xrange(1, int(sqrt(max_n)) / 2 + 1)):
p = k * 2 + 1
for n in xrange(k * 3 + 1, L, p):
a[n] = False
return [ 2 ] + map(lambda k: k * 2 + 1,
ifilter(lambda k: a[k], xrange(1, L)))def gen_phi(n, k0 = 0):
yield 1, 1
for k in takewhile(lambda k: primes[k] <= n, xrange(k0, L)):
p = primes[k]
q = p
m = n
while m >= p:
m /= p
for r, phi in gen_phi(m, k + 1):
yield r * q, q * (p - 1) / p * phi
q *= pdef calc_phi(max_n):
phi = [ 0 ] * (max_n + 1)
for n, p in gen_phi(max_n):
phi[n] = p
return phidef S(n):
if n in memo:
return memo[n]
s = n * (n + 1) / 2
max_m = (n - 1) / M
for m in xrange(2, max_m + 1):
s -= S(n / m)
for k in xrange(1, n / (max_m + 1) + 1):
s -= phi[k] * (n / k - max_m)
memo[n] = s
return sN = 10 ** 6
M = int(sqrt(N) * 1.1)
primes = sieve(M)
L = len(primes)
phi = calc_phi(M)
memo = { }
print S(N) - 1
