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