Project Euler 127

http://projecteuler.net/index.php?section=problems&id=127

a = 1の場合は簡単です。その他の場合を考えましょう。
例えばcが16とすると、c / rad(c) = 8となります。だから、rad(a) * rad(b) < 8となります。しかし、aとbはそれぞれ2以外の因子を持たなければなりません。2以外で最も小さい素数は3、5なので、rad(a) * rad(b) ≥ 15となってcが16のときは解がありません。こうやってcを絞ります。
しかし、こうやってやってもまだ遅いです。aかbをc/2くらい走査しないといけないからです。例えば、c = 64のとき、c / rad(c) = 32だから、rad(a) * rad(b) < 32でなければなりません。そうすると因子の選び方が限られます。3と5、3と7だけです。3と5だとすると、rad(a) = 3ならa < c / 2だから、a = 3, 9, 27、という具合に絞れてきます。しかし、これでも10分近くかかりました。

from itertools import *
from fractions import gcd
import time

def iterate(f, x):
    while True:
        yield x
        x = f(x)

def div_pow(n, p):
    e = 0
    while n % p == 0:
        e += 1
        n /= p
    return e, n

def factorize(max_n):
    a = range(max_n + 1)
    b = [ () for n in a ]
    for p in takewhile(lambda p: p * p <= max_n,
                ifilter(lambda n: a[n] == n, xrange(2, max_n + 1))):
        for n in xrange(p * 2, max_n + 1, p):
            e, a[n] = div_pow(a[n], p)
            b[n] = b[n] + ((p, e),)
    
    for n in ifilter(lambda n: a[n] > 1, xrange(2, max_n + 1)):
        b[n] = b[n] + ((a[n], 1),)
    
    return b

def gen_pairs(a):
    if len(a) == 0:
        yield (), ()
    else:
        e = a[0]
        for p, q in gen_pairs(a[1:]):
#           yield p, q
            yield (e,) + p, q
            yield p, (e,) + q

def rad(n):
    return reduce(lambda x, y: x * y[0], fs[n], 1)

def gen_prime_index_less_c(c, k0):
    return ifilter(lambda k: c % primes[k] != 0, count(k0))

def gen_ps(c):
    def gen_primes_limit(c, n, k = 0):
        for l in takewhile(lambda l: primes[l] <= n,
                                gen_prime_index_less_c(c, k)):
            p = primes[l]
            yield (p,)
            for ps in gen_primes_limit(c, n / p, l + 1):
                yield (p,) + ps
    
    n = c / rad(c)
    return ((x, y) for ps in gen_primes_limit(c, n)
                    for x, y in gen_pairs(ps) if x != () and y != ())

def gen_c(max_n):
    def gen_a(ps, limit):
        if len(ps) == 0:
            yield 1
        else:
            p = ps[0]
            for n in takewhile(lambda n: n <= limit,
                            iterate(lambda x: x * p, p)):
                for a in gen_a(ps[1:], limit / n):
                    yield n * a
    
    for c in xrange(2, max_n + 1):
        for ps_a, ps_b in gen_ps(c):
            for a in gen_a(ps_a, (c - 1) / 2 + 1):
                for b in gen_a(ps_b, c - a):
                    if a + b == c:
                        yield c

def gen_c1(max_n):
    for c in xrange(2, max_n + 1):
        b = c - 1
        if rad(b) * rad(c) < c:
            yield c

t = time.clock()
N = 120000
fs = factorize(N)
primes = filter(lambda n: len(fs[n]) == 1 and fs[n][0][1] == 1,
                                                xrange(2, N + 1))
print sum(gen_c(N)) + sum(gen_c1(N))
print time.clock() - t