Project Euler 153(2)

素因数分解ができればガウス整数での約数の和が求められるようにしましょう。
まず、素数のべき乗について考えます。例えば、3はガウス整数でも素数なので9なら約数は(1, 3, 9)です。一般に4の剰余が3なら同じようになります。
5 = (2 + i)(2 - i)なので、5eの約数は、

(2 + i)e1(2 - i)e2 (0 ≤ e1, e2e)

となります。ここで出てくる約数は必ずしも実部が正になりませんが、和を取るときに90度単位で回転して第1象限または正の整数にします。最後に、2のべきは、2 = (1 + i)(1 - i)ですが、これらは90度ずれているだけなので同じ因数とみなせます。よって、2eの約数は

(1 + i)e1 (0 ≤ e1e)

残りの合成数は、これらの項をもれなく掛け合わせるだけです。
和は、回転して正の整数になるならそれを、そうでなければ、例えば2 + iなら-90度回転させた1 - 2iも実部が正なので、結局実部と虚部の絶対値の和になります。
しかし、この方法では105でもとても遅いです。

from itertools import *
from math import sqrt

def sieve(N):
    a = range(N + 1)
    for p in takewhile(lambda p: p * p <= N,
                    (k for k in xrange(2, N + 1) if a[k] == k)):
        for k in xrange(p * 2, N + 1, p):
            if a[k] == k:
                a[k] = p
    return a

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

def factorize(n):
    if n == 1:
        return ()
    p = fac[n]
    e, n = div_pow(n, p)
    return ((p, e),) + factorize(n)

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 sum_divisors(n):
    def mul(divs1, divs2):
        return [ z1 * z2 for z2 in divs2 for z1 in divs1 ]
    
    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))
    
    a = reduce(lambda x, (p, e): mul(x, prime_pow_divisors(p, e)),
                                                factorize(n), (1,))
    return sum(s(z) for z in a)

def gen_divisors(N):
    for m in xrange(1, N + 1):
        yield 1

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

N = 10 ** 5
fac = sieve(N)
fac_g = fac_gauss(N)
print sum(sum_divisors(n) for n in xrange(1, N + 1))