Project Euler 131

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

n3 + n2p = r3

と書けます。
まず、nが因子pe個含むとしましょう。

n = pe m
n3 + p n2 = p2e+1(pe-1m + 1)m2

最初の項以外はpを含まないので2e+1は3で割り切れます。そうすると、e-1も3で割り切れます。そうすると最初の項以外は、

(pk m)3 + m2

これは、(pk m)3と(pk m + 1)3の間にあるので、立方にはなり得ません。すなわち、npを含みません。
npでないある因子qe個含むとします。

n = qe m
n3 + p n2 = q2em2(qe m + p)

qe m + pqで割り切れないから、eは3の倍数でなければなりません。すなわち、nは立方数となります。

n = m3
n3 + p n2 = m6(m3 + p)
m3 + p = r3

と書けるから、

(r - m)(r2 + r m + m2) = p

左辺の2項のどちらかが1でどちらかがpとなりますが、第2項は1になり得ないので、

r - m = 1
3m2 + 3m + 1 = p
(6m + 3)2 = 12p - 3

丹念に調べると、この式が成り立つのとpが6で割って1余り12p-3が平方数であるのとが同値であることが分かります。

from itertools import *
from math import sqrt

def gen_primes():
    a = [ True ] * L
    for p in (n for n in xrange(2, L) if a[n]):
        for k in xrange(p * 2, L, p):
            a[k] = False
    return (n for n in xrange(2, L) if a[n])

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

L = 10 ** 6
print sum(imap(is_square,
        (12 * p - 3 for p in gen_primes() if p % 6 == 1)))