Project Euler 221

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

Q221.
A = pqr, 1/A = 1/p + 1/q + 1/r
となる正の整数Aの15万番目を求めよ。

整理して、

(p + q)(p + r) = p² + 1

ここで、p > 0, q, r < 0。
右辺を素因数分解してもいいが、やはりコストが大きい。Q216と同じような手法で、エラトステネスのふるい的に素因数分解していく。

def numerize(b, n):

    m = 1

    for p, e in b:

        m *= p ** (n % (e + 1))

        n /= (e + 1)

    return m
def gen_divisors(b):

    def mul(x, y): return x * y

    n = reduce(mul, map(lambda t: t[0] ** t[1], b))

    m = reduce(mul, map(lambda t: t[1] + 1, b))

    for i in xrange(m):

        k = numerize(b, i)

        if k * k <= n:

            yield k
def proc_p(c, p, n):

    if n + p > N:

        return

    for m in xrange(n + p, N + 1, p):

        q = c[m][-1][0]

        e = 0

        while q % p == 0:

            q /= p

            e += 1

        c[m][-1] = (p, e)

        if q != 1:

            c[m].append( (q, 1) )
N = 15 * 10 ** 4

ary = [ ]

set_primes = set()

c = [ [ (n * n + 1, 1) ] for n in xrange(N + 1) ]

for n in xrange(1, N + 1):

    p = c[n][-1][0]

    if p in set_primes:

        continue

    proc_p(c, p, n)

    if p != 2:

        proc_p(c, p, p - n)

    set_primes.add(p)
for p in xrange(1, N + 1):

    n = p * p + 1

    for d in gen_divisors(c[p]):

        q = -p - d

        r = -p - n / d

        A = p * q * r

        ary.append(A)

ary.sort()

print ary[N-1]