\\ maximum square of the modulus mx = 27000 \\ Gaussian integers in the first quadrant, ordered by norm and then by imaginary part quad = concat(vector(sqrtint(mx), r, vector(1+sqrtint(mx-r^2), i, r+(i-1)*I))) mycmp(a,b) = { my (d=norm(a)-norm(b)); if (d==0, d=imag(a)-imag(b) ); sign(d) } quad = vecsort(quad, mycmp) unseen = 1 \\ compute next term other(i) = { my (z=quad[i]); quad[i] = 0; while (quad[unseen]==0, unseen++; ); for (q=unseen, #quad, if (quad[q]!=0, if (z==1 || gcd(z, quad[q])!=1, return (q); ); ); ); print ("# the end"); quit; } { i = 1; for (n=1, 10000, my (z=quad[i]); print (n " " imag(z)); i = other(i); ); } quit