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

import Number
{-- n^3 + n^2 p = n^2 (n+p) = cube
p : prime -> gcd n (n+p) = 1 より
n^2 と n+p は共通因数を持たない
従って n = a^3 , n+p = (a+k)^3 とおける
ここで (a+k)^3 = a^3 + 3 a^2 k + 3 a k^2 + k ^3 = n + k(3a^2+3ak+k^2)
よって、 p = k(3a^2+3ak + k^2)であることが分かる
そして、p素数より、k=1
p = 3a^2 + 3a + 1
--}
p131 lim = takeWhile(<lim). filter isPrime. map p $ [1..]
where p x = 3*x*(x+1)+1
main = print.length.p131$10^6