Project Euler 108

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

Q108.
1/x + 1/y = 1/n の解(x, y > 0)の個数が1000を超える最小のn

(x - n)(y - n) = n²となるから、n²の約数を数えればよい。あとから見るように対象となるnは偶数だから、重複を考えて、約数は2001以上あればよい。
nを1から順に見ていっても求まるが、時間がかかる。逆に素因数分解した形から見ていく。
m = Πp_i^e_iなら、mの約数の個数はΠ(e_i + 1)となる。約数の個数に素数の値自体は関係ないから、素数は小さいものから使い、eも降順とすればよい。n²を素因数分解したとき、e+1は3以上だから、3⁷ = 2187 > 2001だから、素数はせいぜい7つあればよい。すなわち、{ 2, 3, 5, 7, 11, 13, 17 }を使う。分割数のときと同じようにギリギリ2001以上になるeのリストを発生させ、最も小さいnを求める。ただし、使う素数の個数が少ないときに、元の数が非常に大きくなる可能性があるので、そこには配慮して、なるべく大きな数の計算が行われないように工夫した。

from math import log
def make_primes(n):

    primes = [ 2 ]

    m = 3

    while len(primes) < n:

        for p in primes:

            if p * p > m:

                primes.append(m)

                break

            elif m % p == 0:

                break

        m += 2

    return primes
def div(a, b):

    c = { }

    for p, e in a.iteritems():

        c[p] = e

        if b.has_key(p):

            c[p] -= b[p]

    for p, e in b.iteritems():

        if not c.has_key(p):

            c[p] = -e

    return c
def key_with_max_value(a):

    return reduce(lambda x, y: max(x, y, key=lambda x: abs(x[1])),

    a.iteritems())[0]
def less(a, b):

    c = div(a, b)

    pmax = key_with_max_value(c)

    num = 1

    den = 1

    for p, e in c.iteritems():

        if p != pmax:

            if e > 0:

                num *= p ** e

            else:

                den *= p ** -e

    emax = c[pmax]

    if emax > 0:

        while emax > 0:

            num *= pmax

            if num > den:

                return False

            emax -= 1

        return True

    else:

        while emax < 0:

            den *= pmax

            if num < den:

                return True

            emax += 1

        return False
def Min(a, b):

    if less(a, b):

        return a

    else:

        return b
def gen_div(n, m, r = 0):

    if r == 0:

        r = n

    elif r > n:

        r = n

    

    if m == 1:

        if n <= r and n & 1 and n >= 3:

            yield [ n ]

    elif n >= 3:

        for k in xrange(3, r + 1, 2):

            for a in gen_div((n + k - 1) / k, m - 1, k):

                yield [ k ] + a
def exp_to_num(a):

    m = { }

    for i in range(len(a)):

        m[primes[i]] = a[i] / 2

    return m
def min_number_core(M, m):

    return reduce(Min, map(exp_to_num, gen_div(M, m)))
def min_number(M):

    a = min_number_core(M, 1)

    for i in range(2, 8):

        b = min_number_core(M, i)

        a = Min(a, b)

    return reduce(Min, map(lambda i: min_number_core(M, i),

    range(1, max_num_primes + 1)))
def val(a):

    def pow(x): return x[0] ** x[1]

    return reduce(lambda x, y: x * y, map(pow, a.iteritems()))
N = 1000

M = 2 * N + 1

max_num_primes = int(log(M) / log(3)) + 1

primes = make_primes(max_num_primes)
print val(min_number(M))