Project Euler 72(1)

Problem 72
n/dという分数を考える。ここで、ndは正の整数。n<dでHCF(n,d)=1のとき、既約真分数と呼ぶ。
d≤8の既約真分数の集合を昇順に並べると次のようになる:
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
この集合は21個の要素がある。
d≤1000000の既約真分数の集合の要素数は?
http://projecteuler.net/index.php?section=problems&id=72

これは要するに、N=1000000として

φ(2) + φ(3) + … + φ(N)

を求めよ、ということです。φ(n)はオイラーのφ関数です。
さて、φ(n)は求めるのに通常は素因数分解が要ります。

n = p1e1...pmem
φ(n) = p1e1-1(p1 - 1)...pmem-1(pm - 1)

素因数分解を求めるには2つの方法があります。一つはエラトステネスのふるい的方法です(Code1)。


もう一つは次々に素因数分解を生成していくやり方です。N=10で考えましょう。10以下の素数が必要になります。まず、1を返します。次に、素因数分解したときの最小の素数を選びます。2なら、10以下で21, 22,23となります。21を選んだら10/2=5以下で3以上の素数のみを使った数を求めます。最初に1を返して、次に31と51を返します。これで、21,2131,2151という素因数分解が得られました。このようにすると再帰的にN以下の全ての素因数分解が得られます。
今回は素因数分解ではなく直接φを計算して生成しています(Code2)。
今回の問題ではふるい的方法よりこちらのほうが速くなりました。


よく考えると、いちいち個々のφを生成する必要はありません。

φ(2) + φ(6) + φ(10) = φ(2)(1 + φ(3) + φ(5))

右辺第二項は3以上の素数を使った5以下のφの和と言うことができます。これでさらに速くなりました(Code3)。




# Code1
from itertools import ifilter, imap
from math import sqrt

def sieve(max_n):
L = (max_n + 1) / 2
a = range(1, max_n + 1, 2)
for k in ifilter(lambda k: a[k], xrange(1, int(sqrt(max_n)) / 2 + 1)):
p = a[k]
for l in xrange(k + p, L, p):
if a[l] == l * 2 + 1:
a[l] = p
return a

def extract_pow(n, p):
q = p
while n % q == 0:
q *= p
return q / p

def phi(n):
if n == 1:
return 1
elif b[n]:
return b[n]
else:
if n % 2 == 0:
q = extract_pow(n, 2)
m = q / 2 * phi(n / q)
else:
p = a[(n-1)/2]
q = extract_pow(n, p)
m = q * (p - 1) / p * phi(n / q)
b[n] = m
return m

N = 10 ** 6
a = sieve(N)
b = [ 0 ] * (N + 1)
print sum(imap(phi, xrange(2, N + 1)))


# Code2
from itertools import ifilter, imap, takewhile
from math import sqrt

def 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
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 phi in gen_phi(m, k + 1):
yield q * (p - 1) / p * phi
q *= p

N = 10 ** 6
primes = sieve(N)
L = len(primes)
print sum(gen_phi(N)) - 1


# Code3
from itertools import ifilter, imap, takewhile
from math import sqrt

def 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 sum_phi(n, k0 = 0):
s = 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
s += q * (p - 1) / p * sum_phi(m, k + 1)
q *= p
return s

N = 10 ** 6
primes = sieve(N)
L = len(primes)
print sum_phi(N) - 1