%I #37 Aug 18 2017 09:41:43
%S 1,4,7,13,16,19,25,28,31,37,43,49,52,61,64,67,73,76,79,91,97,100,103,
%T 109,112,121,124,127,133,139,148,151,157,163,169,172,175,181,193,196,
%U 199,208,211,217,223,229,241,244,247,256,259,268,271,277,283,289,292
%N Numbers of the form 3*(x^2 + xy + y^2 + x + y) + 1 where x and y are integers.
%C Closed under multiplication.
%C Löschian numbers of the form 3*k+1. - _Altug Alkan_, Nov 18 2015
%H Reinhard Zumkeller, <a href="/A202822/b202822.txt">Table of n, a(n) for n = 1..10000</a>
%H Joerg Arndt, <a href="https://arxiv.org/abs/1607.02433">Plane-filling curves on all uniform grids</a>, arXiv preprint arXiv:1607.02433 [math.CO], 2016.
%F A033685(n) != 0 if and only if n is in the set.
%t nf[{i_,j_}]:=3(i^2+i*j+j^2+i+j)+1; Union[nf/@Tuples[Range[-10,10],2]] (* _Harvey P. Dale_, Dec 31 2011 *)
%o (PARI) isA(n) = if(n%3 == 0, 0, 0 != sumdiv( n, d, kronecker( -3, d)))
%o (PARI) x='x+O('x^500); p=eta(x)^3/eta(x^3); for(n=0, 499, if(polcoeff(p, n) != 0 && n%3==1, print1(n, ", "))) \\ _Altug Alkan_, Nov 18 2015
%o (PARI) is(n)=(n%3==1) && #bnfisintnorm(bnfinit(z^2+z+1), n); \\ _Joerg Arndt_, Jan 04 2016
%o (PARI) list(lim)=my(v=List(), y, t); for(x=0, sqrtint(lim\3), my(y=x, t); while((t=x^2+x*y+y^2)<=lim, if((x-y)%3, listput(v, t)); y++)); Set(v) \\ _Charles R Greathouse IV_, Jul 05 2017
%o (Haskell)
%o a202822 n = a202822_list !! (n-1)
%o a202822_list = filter ((== 1) . flip mod 3) a003136_list
%o -- _Reinhard Zumkeller_, Nov 16 2015
%Y Cf. A033685, A045833.
%Y Subsequence of A003136, A260682 (subsequence).
%K nonn
%O 1,2
%A _Michael Somos_, Dec 25 2011