Project Euler 153（3）

素因数分解をするのではなく、素数の積を生成するのと同じ方法で約数の組を作っていきます。最も素直に書いてみました。
しかし、やっぱり遅いですね。

from itertools import *
from math import sqrt
import time

def sieve(max_n):
    L = max_n / 1
    a = [ True ] * L
    for p in takewhile(lambda n: n * n < L,
                    (n for n in xrange(2, L) if a[n])):
        for k in xrange(p * 2, L, p):
            a[k] = False
    primes = [ n for n in xrange(2, L) if a[n] ]
    
    for offs in xrange(L, max_n + 1, L):
        a = [ True ] * L
        for p in takewhile(lambda n: n * n < offs + L, primes):
            for m in xrange((offs + p - 1) / p * p, offs + L, p):
                a[m-offs] = False
        primes.extend(k + offs for k in xrange(L)
                                if a[k] and k + offs <= max_n)
    return primes

def prime_pow_divisors(p, e):
    if p % 4 == 3:
        return (p ** e1 for e1 in xrange(e + 1))
    elif p % 4 == 1:
        z1 = fac_g[p]
        z2 = z1.conjugate()
        return (z1 ** e1 * z2 ** e2 for e1 in xrange(e + 1)
                                    for e2 in xrange(e + 1))
    else:
        return ((1 + 1j) ** e1 for e1 in xrange(e * 2 + 1))

def mul_divs(divs1, divs2):
    return [ z1 * z2 for z1 in divs1 for z2 in divs2 ]

def gen_divisors(n, k0 = 0):
    yield (1,)
    if k0 < len(primes) and n >= primes[k0]:
        g_primes = takewhile(lambda p: p <= n, islice(primes, k0, None))
        for k, p in enumerate(g_primes, k0):
            for e in takewhile(lambda e: p ** e <= n, count(1)):
                divs1 = list(prime_pow_divisors(p, e))
                for divs2 in gen_divisors(n / p ** e, k + 1):
                    yield mul_divs(divs1, divs2)

def sum_divisors(divs):
    def s(z):
        if z.imag == 0:
            return int(abs(z.real))
        elif z.real == 0:
            return int(abs(z.imag))
        else:
            return abs(int(z.real)) + abs(int(z.imag))
    
    return sum(s(z) for z in divs)

def fac_gauss(N):
    M = int(sqrt(N))
    fac_g = [ 0 ] * (N + 1)
    for a in xrange(2, M + 1):
        for b in xrange(1, a):
            n = a * a + b * b
            if n <= N:
                fac_g[n] = a + b * 1j
    return fac_g

t = time.clock()
N = 10 ** 5
primes = sieve(N)
fac_g = fac_gauss(N)
print sum(imap(sum_divisors, gen_divisors(N)))
print time.clock() - t