http://projecteuler.net/index.php?section=problems&id=131
n3 + n2p = r3
と書けます。
まず、nが因子pをe個含むとしましょう。
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の間にあるので、立方にはなり得ません。すなわち、nはpを含みません。
nはpでないある因子qをe個含むとします。
n = qe m
n3 + p n2 = q2em2(qe m + p)
qe m + pはqで割り切れないから、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)))