Project Euler 153（4）

例えば、5の実部が正の約数は、

1, 2 + i, 1 - 2i, 2 - i, 1 + 2i, 5

ですが、15の約数は、5の約数とそれを3倍したものになります。すなわち、s(15) = 4s(15)となります。105なら、s(105) = (1 + 3 + 7 + 21)s(5)となります。
ガウス整数の約数を計算しなければならないのは、4の剰余が3でない素数の積のみで、その和を取ったあとに4の剰余が3の素数の積の約数の和をかければよいです。こうすることにより高速化を図れます。

from itertools import *
from math import sqrt

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 gen_pow(n, e):
    m = 1
    yield 1
    for k in xrange(e):
        m *= n
        yield m

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

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

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

def sum_divisors((n, 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) * sum_divisors3(N / n)

def sum_divisors3(n, k0 = 0):
    s = 1
    if k0 < len(primes) and n >= primes[k0]:
        g_primes = takewhile(lambda p: p <= n, islice(primes, k0, None))
        for k, p in ((k, p) for k, p in enumerate(g_primes, k0) if p % 4 == 3):
            for e in takewhile(lambda e: p ** e <= n, count(1)):
                m = p ** e
                s += (m * p - 1) / (p - 1) * sum_divisors3(n / m, k + 1)
    return s

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
primes = sieve(N)
fac_g = fac_gauss(N)
print sum(imap(sum_divisors, gen_divisors1(N)))