isgaussianprime(z) = { my (x); if (real(z)==0, ((x=abs(imag(z)))%4)==3 && isprime(x), imag(z)==0, ((x=abs(real(z)))%4)==3 && isprime(x), isprime(norm(z))); } { T = matrix(11,11); n=0; z=0; dir=1; for (w=1, oo, for (t=0, 1, for (s=1, w, n++; if (isgaussianprime(z), print (k++ " " n); if (k==10 000, break (3)); ); z+=dir; ); dir*=I; ); ); } quit