Project Euler 141

http://projecteuler.net/index.php?section=problems&id=141

rがdとqの間のとき、公差をaとして、

n = r + r/a * ra = r + r²

これは平方数になり得ないので、rが一番小さいということになります。そうすると、

n = r + ra * ra² = r + r²a³

aをs/t（sとtは互いに素）として、

n = r + r²s³ / t³

r² / t³が整数だから、r = kt²と表せて、

n = kt² + k²s³t

これを素直にコーディングしても40秒以内に答えが出ました。
しかし、例えばt = 1とすると、

n = m² = k + k²s³
4s⁶k + 4s³k = 4s³m²
(2s³k + 1) - 4s(s m)² = 1

と、ペル方程式になります。
しかし、初期解を求めるのに時間がかかり、かえって遅くなってしまいました。この件は保留にしておきます。

from itertools import *
from fractions import gcd
from math import sqrt

def is_square(n):
    m = int(sqrt(n))
    return m * m == n

def gen_progressive_perfect_squares(N):
    for s in takewhile(lambda s: 1 + s * s * s < N, count(2)):
        s3 = s * s * s
        for t in ifilter(lambda t: gcd(s, t) == 1,
                    takewhile(lambda t: t * (t + s3) < N, xrange(1, s))):
            for n in ifilter(is_square, takewhile(lambda n: n < N,
                            (k * t * (t + s3 * k) for k in count(1)))):
                yield n

N = 10 ** 6
print sum(gen_progressive_perfect_squares(N))