Project Euler 75

Problem 75
（前略）Lをワイヤーの長さとして、L ≤ 1,500,000が唯一つの辺の長さが整数の直角三角形を成すLの個数を求めよ。
http://projecteuler.net/index.php?section=problems&id=75

ピタゴラス数の生成式を使うと周の長さは

a + b + c = 2lm(m+n)

となり、m + n -> nと置き換えて、

L = 2lmn
m ≤ n ≤ 2m
m, nは互いに素
nは奇数

となります。だから、Lを素因数分解して約数を出して条件に合うものを選べば、Lになる直角三角形を数えることができます（Code1）。ただし、非常に遅いです。
素因数分解するのではなく、条件に合うl,m,nからLを出して個数をカウントしていくほうが速いです。lを最後にするのが速いようです（Code2）。

# Code1
from itertools import ifilter, imap
from fractions import gcd

def sieve(max_n):
    a = range(max_n + 1)
    for p in ifilter(lambda n: a[n] == n, xrange(2, max_n + 1)):
        for n in ifilter(lambda n: a[n] == n, xrange(p * 2, max_n + 1, p)):
            a[n] = p
    return a

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

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

def gen_divisors(f, k = 0):
    if k == len(f):
        yield 1, 1
    else:
        p, e0 = f[k]
        for d1, d2 in gen_divisors(f, k + 1):
            for e in range(e0 + 1):
                yield d1 * p ** e, d2 * p ** (e0 - e)

def gen_divisors2(f, k = 0):
    if k == len(f):
        yield 1, []
    else:
        p, e0 = f[k]
        for d1, f2 in gen_divisors2(f, k + 1):
            for e in range(e0 + 1):
                yield d1 * p ** e, [(p, e0 - e)] + f2

def gen_divs(f):
    for d1, f2 in gen_divisors2(f):
        m = 1
        if len(f2) > 0 and f2[0][0] == 2:
            m = 1 << f2[0][1]
            f2 = f2[1:]
        for d2, d3 in gen_divisors(f2):
            yield d1, m * d2, d3

def is_valid(d):
    l, m, n = d
    return m < n < m * 2 and gcd(m, n) == 1

def count_right_angle_triangles(n):
    return sum(imap(is_valid, gen_divs(factorize(n / 2))))

L = 15 * 10 ** 5
M = L / 2
a = sieve(M)
print sum(imap(lambda n: n == 1, imap(
            count_right_angle_triangles, xrange(2, L + 1, 2))))

# Code2
from itertools import imap
from fractions import gcd
from math import sqrt

t = time.clock()
L = 15 * 10 ** 5
M = L / 2
a = [ 0 ] * (M + 1)
for m in xrange(2, int(sqrt(M)) + 1):
    for n in xrange(m + 2 if m & 1 else m + 1, min(M / m + 1, m * 2), 2):
        if gcd(m, n) == 1:
            p = m * n
            for r in xrange(p, M, p):
                a[r] += 1

print sum(imap(lambda n: a[n] == 1, xrange(M + 1)))