Project Euler 248

プロジェクトオイラー
http://projecteuler.net/

Q248.
φ(n) = 13!となる小さいほうから15万番目のn。

n = p^eなら、φ(n) = (p - 1)p^e-1で、かつ乗法的である（(x, y) = 1なら、φ(xy) = φ(x)φ(y)）。
まず、1を引くと13!の約数になる素数を全て求める。その組合せで、再帰的に13!になる整数を全て求め、ソートして15万番目の整数を表示する。もっと細かく制御すれば速くなるんだろうけど。

from itertools import count

from math import sqrt, factorial
def is_prime(n):

    for p in primes:

        if p * p > n:

            return True

        elif n % p == 0:

            return False

    return True
def make_primes(n):

    primes.append(2)

    m = int(sqrt(n + 0.5))

    for p in xrange(3, m + 1, 2):

        if is_prime(p):

            primes.append(p)

    

    a = [ 0 ] * ((n + 15) / 16)

    for p in primes:

        if p * p > n:

            break

        if p == 2:

            continue

        for k in xrange(3 * p, n, 2 * p):

            a[k>>4] |= 1 << (k & 15)

    

    for k in xrange((m + 1) / 2 * 2 + 1, n + 1, 2):

        if (a[k>>4] & (1 << (k & 15))) == 0:

            primes.append(k)
def factorize(n):

    s = 0

    for k in count(0):

        p = primes[k]

        if p > N:

            break

        bit = 1 << k * 4

        while n % p == 0:

            n /= p

            s += bit

    if n == 1:

        return s

    else:

        return 0
def gen_factors(f, k = 0):

    p = primes[k]

    e_max = (f >> 4 * k) & 15

    if f >> 4 * (k + 1) == 0:

        for e in range(e_max + 1):

            yield p ** e

    else:

        for n in gen_factors(f, k + 1):

            for e in range(e_max + 1):

                yield p ** e * n
def make_factors(f):

    a = [ ]

    for n in gen_factors(f):

        m = n + 1

        if is_prime(m):

            a.append( (m, factorize(n)) )

    return a
def cmp(f1, f2):

    b = 0

    for n in xrange(6):

        e1 = (f1 >> n * 4) & 15

        e2 = (f2 >> n * 4) & 15

        if e1 > e2:

            return -1

        elif e1 < e2:

            b = 1

    return b
def gen_number(factors, f0, k0 = 0):

    if k0 < len(factors):

        for k in xrange(k0, len(factors)):

            p = factors[k][0]

            q = p

            f = f0

            f1 = factors[k][1]

            if p <= N:

                fp = factorize(p)

            while True:

                c = cmp(f1, f)

                if c == 0:

                    yield q

                    if p == 2 and q == p:

                        yield p * 2

                    break

                elif c == 1:

                    f -= f1

                    for n in gen_number(factors, f, k + 1):

                        yield q * n

                    if p > N:

                        break

                    q *= p

                    f1 = fp

                else:

                    break
def make_number_array(f):

    a = [ ]

    factors = make_factors(f)

    for n in gen_number(factors, f):

        a.append(n)

    return a
N = 13

L = factorial(N)

M = 15 * 10 ** 4

primes = [ ]

make_primes(int(sqrt(L + 1)))

F = factorize(L)

ary = make_number_array(F)

ary.sort()

print ary[M-1]