ガウス整数

ガウス整数は、a, bを整数としたとき、a + ibと表せる数を言います。
ガウス整数は、素因数分解の一意性があります。例えば5なら、

5 = (2 + i)(2 - i)

と一意に分解できます。ただし、-1、±iを掛けた数は同じ因数とみなします。すなわち、

5 = (-1 + 2i)(-1 - 2i)
5 = (-2 - i)(-2 + i)
5 = (1 - 2i)(1 + 2i)

も同じ分解とみなします。
この性質は、原点を中心とする円周上の格子点の個数を数えるときに使えます。例えば円周を

x² + y² = 65

とします。

1² + 2² = (1 + 2i)(1 - 2i) = 5
2² + 3² = (2 + 3i)(2 - 3i) = 13

なので、

(1 + 2i)(2 + 3i) = -4 + 7i
(1 + 2i)(2 - 3i) = 8 + i
(1 - 2i)(2 + 3i) = 8 - i
(1 - 2i)(2 - 3i) = -4 - 7i

の4つが円周上の格子点になります。例えば上2つで確かめると、

(-4)² + 7² = 16 + 49 = 65
8² + 1² = 64 + 1 = 65

しかし、i, -1, -iを掛けた座標も90度回っただけなので、全部で16個の格子点があることになります。1象限につき4つになります。第1象限では、

(4, 7), (7, 4), (8, 1), (1, 8)

になります。
一般に素数pが4m + 1と表せるときにx² + y² = pとなるx, yがあります。ガウス整数の素因数分解の一意性のためこの解は対称の位置を除いて一つです。素数pが4m + 3と表されるときに解がないことは明らかです。
よって、nが異なるk個の素数の積で表されるとき、それが全て4m + 1型なら、第1象限に2^k個あることになります。一つでも4m + 3型の素数が含まれるなら、解はありません。
nがp^eと表される場合は、共役のガウス整数を幾つずつ取るかということになるので、e + 1個となります。例えば125なら、

(1 + 2i)³ = -11 - 2i
(1 + 2i)²(2 + i) = -10 + 5i
(1 + 2i)(2 + i)² = (1 + 2i)(3 + 4i) = -5 + 10i
(2 + i)³ = 2 + 11i

の4つになります。また、pが4m + 3型でeが偶数なら一つの解があります。
さて、最後にnが2を含む場合です。

2 = (1 + i)(1 - i)

ですが、

(1 - i)i = 1 + i

なので、すなわち90度回転分違うだけなので、2つは同じ因子だとみなされます。すなわち、2が含まれていても含まれていなくても解の個数に影響はないことになります。
まとめると、

n = p₁^e₁ … p_m^e_m

ならx² + y² = n上の第1象限の格子点の数は、

g(p₁, e₁) … g(p_m, e_m)

ここでgは、

g(p, e) = 1 (p = 2)
g(p, e) = e + 1 (p = 4k + 1)
g(p, e) = 1 (p = 4k + 3, e : even)
g(p, e) = 0 (p = 4k + 3, e : odd)

ただし、第1象限という意味では軸上は含まないので、nが平方数の場合は1小さくなります。

from itertools import *

def factorize(n):
    def div_pow(n, p):
        e = 0
        while n % p == 0:
            e += 1
            n /= p
        return e, n
    
    f = ()
    for k in takewhile(lambda k: k * k <= n, count(2)):
        e, n = div_pow(n, k)
        if e > 0:
            f = f + ((k, e),)
    if n > 1:
        f = f + ((n, 1),)
    return f

def g((p, e)):
    if p == 2:
        return 1
    elif p % 4 == 1:
        return e + 1
    else:
        return 1 if e % 2 == 0 else 0

def is_square(f):
    return all(e % 2 == 0 for (p, e) in f)

def num_lattice_points(n):
    f = factorize(n)
    return reduce(lambda x, y: x * g(y), f, 1) \
                        - (1 if is_square(f) else 0)

for n in xrange(10000, 10030):
    print n, num_lattice_points(n)

Project Eulerに必要な用語集
http://d.hatena.ne.jp/inamori/20091225/p1