primosConCubos :: [Integer]

λ> take 6 primosConCubos
[7,19,37,61,127,271]
λ> length (takeWhile (< 1000000) primosConCubos)
173

import Data.Numbers.Primes

-- 1ª solución
primosConCubos :: [Integer]
primosConCubos = [p | x <- [1..],
                      n <- [1..x],
                      (x^3 - n^3) `mod` (n^2) == 0,
                      let p = (x^3 - n^3) `div` (n^2),
                      isPrime p]

-- Análisis del problema
-- =====================

-- El número p es primo con cubo si existen dos enteros positivos n y k
-- tales
--    n^3 + n^2p = k^3
-- Se tienen las siguientes equivalencias
--    n^3 + n^2p = k^3
--    <=> n^3(1 + p/n) = k^3
--    <=> n^3((n+p)/n) = k^3
--    <=> n(((n+p)/n)^(1/3)) = k
--    <=> n ((n+p)^(1/3)/(n^(1/3))) = k
-- Por tanto, n y n+p son cubos perfectos; es decir, existen enteros
-- positivos x e y tales que
--    n = x^3                               (1)
--    n+p = y^3                             (2)
-- De (2),
--    p = y^3 - n                           (3)
-- Sustituyendo (1) en (3),
--    p = y^3 - x^3                         (4)
-- Simplificando
--    p = (y-x)(x^2 + xy + y^2)
-- de donde, por ser p primo, se tiene que
--    y-x = 1                               (5)
--    p = x^2 + xy + y^2                    (6)
-- De (4) y (5), se tiene que p es la diferencia de dos cubos
-- consecutivos. Además, despejando y en (5) y sustituyendo en (6), se
-- tiene que
--    p = x^2 + x(1+x) + (1+x)^2
--      = x^2 + x + x^2 + 1 + 2x + x^2
--      = 3x^2+ 3x + 1

-- 2ª solución
-- ===========

-- Como primos cubanos

-- 3ª solución
-- ===========

primosConCubos3 :: [Integer]
primosConCubos3 = [p | x <- [1..],
                       let p = 3*x^2+ 3*x + 1,
                       isPrime p]