Project Euler 110

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

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

108とほとんど同じ。少し工夫した。
例えば、素因数分解したときに素数が2つというのはありえない。約数の個数が100として、2か3のどちらかの指数が9以上となる。より厳しい2の指数が9以上とすると、同じ約数の個数でも、2⁹ > 2⁴5¹だから、素数が3つ以上あったほうが小さくなる。このような考えに基づくと、素数は8個以上必要であることがわかる。

from itertools import count

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):

    if a == 0:

        return False

    elif b == 0:

        return True

    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):

    print M, m

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

    min_num_primes = calc_min_num_primes(M)

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

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

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

    return reduce(lambda x, y: x * y, map(pow, a.iteritems()))
def calc_min_num_primes(M):

    for n in count(2):

        min_e = int(M ** (1.0 / n))

        if (min_e + 1) / 2 * log(2) < log(primes[n]):

            return n
N = 4 * 10 ** 6

M = 2 * N + 1

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

primes = make_primes(max_num_primes)

print max_num_primes
print val(min_number(M))