OFFSET
1,1
COMMENTS
Equivalently, numbers of the form (x+y)*(x^2 + y^2) where x and y are coprime positive integers.
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
EXAMPLE
For x=1, y=1, x^3+x^2*y+x*y^2+y^3 = 4, so 4 is in the sequence.
For x=1, y=2, x^3+x^2*y+x*y^2+y^3 = 15, so 15 is in the sequence.
For x=2, y=2, x^3+x^2*y+x*y^2+y^3 = 32, but GCD(x,y)<>1, so 32 is not in the sequence.
MAPLE
N:= 10000: # for terms <= N
S:= {}:
for x from 1 while (x+1)*(x^2+1) < N do
V:= select(`<=`, map(y -> (x+y)*(x^2+y^2), select(y -> igcd(x, y)=1, {seq(i, i=1..min(x, (N-x^3)/x^2))})), N);
S:= S union V;
od:
sort(convert(S, list)); # Robert Israel, Sep 21 2020
MATHEMATICA
max = 5000; T0 = {}; xm = Ceiling[Sqrt[max]];
While[T = T0;
T0 = Table[x^3 + x^2 y + x y^2 + y^3, {x, 1, xm}, {y, x, xm}] //
Flatten // Union // Select[#, # <= max &] &; T != T0, xm = 2 xm];
T (* T=A336607 *)
(* Now, exclude a(n) such that a(n)=k^3*a(m) for m<n and k>=2 is integer *)
T2 = T; n = 1;
While[n <= Length[T2],
t1 = T2[[n]]; t2 = Last[T2]; max2 = 1 + (t2/t1)^(1/3);
T2 = Complement[T2, Table[t1*k^3, {k, 2, max2}]];
n++;
];
T2 (* T2=A336995 *)
PROG
(PARI) upto(limit)={my(L=List(), b=sqrtnint(limit, 3)); for(x=1, b, for(y=1, b, my(t=(x+y)*(x^2+y^2)); if(t<=limit && gcd(x, y)==1, listput(L, t)) )); Set(L)}
upto(4000) \\ Andrew Howroyd, Aug 10 2020
CROSSREFS
KEYWORD
nonn,easy
AUTHOR
César Eliud Lozada, Aug 10 2020
STATUS
approved