cp's OEIS Frontend

A336995 Numbers of the form x^3 + x^2*y + x*y^2 + y^3, where x and y are coprime positive integers.

%I A336995 #21 Jan 11 2025 13:47:45
%S A336995 4,15,40,65,85,156,175,203,259,272,369,400,477,580,585,671,715,803,
%T A336995 820,888,935,1105,1111,1157,1261,1417,1464,1484,1625,1695,1820,1885,
%U A336995 2055,2080,2336,2380,2465,2533,2595,2669,2848,2873,2955,3060,3145,3439,3485,3492
%N A336995 Numbers of the form x^3 + x^2*y + x*y^2 + y^3, where x and y are coprime positive integers.
%C A336995 Equivalently, numbers of the form (x+y)*(x^2 + y^2) where x and y are coprime positive integers.
%H A336995 Robert Israel, <a href="/A336995/b336995.txt">Table of n, a(n) for n = 1..10000</a>
%e A336995 For x=1, y=1, x^3+x^2*y+x*y^2+y^3 = 4, so 4 is in the sequence.
%e A336995 For x=1, y=2, x^3+x^2*y+x*y^2+y^3 = 15, so 15 is in the sequence.
%e A336995 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.
%p A336995 N:= 10000: # for terms <= N
%p A336995 S:= {}:
%p A336995 for x from 1 while (x+1)*(x^2+1) < N do
%p A336995    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);
%p A336995    S:= S union V;
%p A336995 od:
%p A336995 sort(convert(S,list)); # _Robert Israel_, Sep 21 2020
%t A336995 max = 5000; T0 = {}; xm = Ceiling[Sqrt[max]];
%t A336995 While[T = T0;
%t A336995 T0 = Table[x^3 + x^2 y + x y^2 + y^3, {x, 1, xm}, {y, x, xm}] //
%t A336995      Flatten // Union // Select[#, # <= max &] &; T != T0, xm = 2 xm];
%t A336995 T (* T=A336607 *)
%t A336995 (* Now, exclude a(n) such that a(n)=k^3*a(m) for m<n and k>=2 is an integer *)
%t A336995 T2 = T; n = 1;
%t A336995 While[n <= Length[T2],
%t A336995   t1 = T2[[n]]; t2 = Last[T2]; max2 = 1 + (t2/t1)^(1/3);
%t A336995   T2 = Complement[T2, Table[t1*k^3, {k, 2, max2}]];
%t A336995   n++;
%t A336995   ];
%t A336995 T2 (* T2=A336995 *)
%o A336995 (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)}
%o A336995 upto(4000) \\ _Andrew Howroyd_, Aug 10 2020
%Y A336995 Subsequence of A336607.
%K A336995 nonn,easy
%O A336995 1,1
%A A336995 _César Eliud Lozada_, Aug 10 2020