Project Euler 72（1） - inamori’s diary

Problem 72
n/dという分数を考える。ここで、nとdは正の整数。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 = p₁^e₁...p_m^e_m
φ(n) = p₁^e₁-1(p₁ - 1)...p_m^e_m-1(p_m - 1)

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

もう一つは次々に素因数分解を生成していくやり方です。N=10で考えましょう。10以下の素数が必要になります。まず、1を返します。次に、素因数分解したときの最小の素数を選びます。2なら、10以下で2¹, 2²,2³となります。2¹を選んだら10/2=5以下で3以上の素数のみを使った数を求めます。最初に1を返して、次に3¹と5¹を返します。これで、2¹,2¹3¹,2¹5¹という素因数分解が得られました。このようにすると再帰的に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