Project Euler 108

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

(x - n)(y - n) = n²

だから、nを素因数分解して、

n = p₁^e₁…p_m^e_m

とすると、解の個数は、

((2e₁ + 1)…(2e_m + 1) + 1) / 2

となります。
単に割っていって素因数分解するのは遅いので、エラトステネスのふるい的に素因数分解をします。どこまで素因数分解すればいいのかわからないので、少しずつまとめて計算します。ベンチマークしてみると、500ずつが最も速いという結果になりました。
本当は素因数分解をする必要も無いのですが、それは次で。

from itertools import *

def sieve(max_n):
    a = [ True ] * (max_n + 1)
    for p in (n for n in xrange(2, max_n + 1) if a[n]):
        for k in xrange(p * 2, max_n + 1, p):
            a[k] = False
    return [ n for n in xrange(2, max_n + 1) if a[n] ]

def gen_factorize():
    def div_pow(n, p):
        e = 0
        while n % p == 0:
            e += 1
            n /= p
        return e, n
    
    start = 2
    while True:
        a = range(start, start + L)
        b = [ [] for k in xrange(start, start + L) ]
        for p in takewhile(lambda p: p * p < start + L, primes):
            if p >= start:
                begin = p * 2
            else:
                begin = (start + p - 1) / p * p
            for n in xrange(begin, start + L, p):
                k = n - start
                e, a[k] = div_pow(a[k], p)
                b[k].append((p, e))
        
        for k in xrange(L):
            n = start + k
            if len(b[k]) == 0:
                b[k] = [ (n, 1) ]
                primes.append(n)
            elif a[k] > 1:
                b[k].append((a[k], 1))
        
        for f in b:
            yield f
        
        start += L

def count_solutions(f):
    def c(f, k):
        if k == len(f):
            return 1
        else:
            return (f[k][1] * 2 + 1) * c(f, k + 1)
    
    return (c(f, 0) + 1) / 2

N = 1000
L = 500
primes = sieve(L)
print dropwhile(lambda x: x[1] <= N,
    izip(count(2), imap(count_solutions, gen_factorize()))).next()[0]