list(lim)= { my(u=powers(2,logint(lim\=1,2)), v=List(), p1=v, p3=v, pr, t=1); \\ Products of primes = 1 mod 4 forprime(p=5,, if(p%4>1, next); t*=p; if(t>lim, break); listput(p1,t) ); \\ Products of primes = 3 mod 4 t=1; forprime(p=3,, if(p%4<3, next); t*=p; if(t>lim, break); listput(p3,t) ); \\ Append products of primes = 1 mod 4 for(i=1,#p1, pr=1; for(e=1,logint(lim,p1[i]), pr*=p1[i]; for(j=1,#u, t=pr*u[j]; if(t>lim, break); listput(v,t) ) ); if(p1[i]^2lim, break); listput(v,t) ) ); if(p3[i]^2