%I #20 Nov 24 2020 03:26:08
%S 13,25,37,52,61,73,97,100,109,117,148,157,169,181,193,208,225,229,241,
%T 244,277,289,292,313,325,333,337,349,373,388,397,400,409,421,433,436,
%U 457,468,481,541,549,577,592,601,613,625,628,637,657,661,673,676,709,724,733
%N Löschian numbers (A003136) which are the sum of 2 nonzero squares.
%C n is in the sequence iff 4*n is.
%C If a(n) is a prime number, a(n) mod 12 = 1.
%C Prime terms of sequence are listed in A068228 that lists generalized cuban primes (A007645) which are the sum of 2 nonzero squares.
%C Also positive numbers of the form x^2 - 3*y^2 (A084916) that are the sum of 2 nonzero squares. - _Frank M Jackson_, Oct 13 2019
%H Amiram Eldar, <a href="/A264732/b264732.txt">Table of n, a(n) for n = 1..10000</a>
%e a(1) = 13 because 13 = 3^2 + 3*1 + 1^2 = 3^2 + 2^2.
%e a(2) = 25 because 25 = 5^2 + 5*0 + 0^2 = 4^2 + 3^2.
%e a(3) = 37 because 37 = 4^2 + 4*3 + 3^2 = 6^2 + 1^2.
%t Select[Range@750, Length[PowersRepresentations[#, 2, 2] /. {0, _}->Nothing]>0 && Reduce[#==x^2+x*y+y^2, {x, y}, Integers]=!=False &] (* _Frank M Jackson_, Oct 13 2019 *)
%o (PARI) isok(n) = { for( i=1, #n=factor(n)~%4, n[1, i]==3 && n[2, i]%2 && return); n && ( vecmin(n[1, ])==1 || (n[1, 1]==2 && n[2, 1]%2))}
%o is(n) = #bnfisintnorm(bnfinit(z^2+z+1), n);
%o for(n=1, 1e3, if( is(n) && isok(n), print1(n, ", ")))
%Y Cf. A000404, A003136, A068228, A084916.
%K nonn
%O 1,1
%A _Altug Alkan_, Nov 22 2015