Can we find Gaussian primes $\pi = 1 + 8 \mathbb{Z}[i]$ with $N(\pi) < 1000$? It's an exercise in computational number theory.  Either by hand or by computer, can we find the Gaussian primes $\pi = 1 + 8 \mathbb{Z}[i]$?  To keep the list finite I guess we could have $N(\pi) < 1000$.
For example, is $\mathfrak{p} = (1+8i)$ a prime?  Or $\mathfrak{p} = (-7 + 8i)$? I don't even know how to index the primes less than these.  The norms are $1^2 + 8^2 = 65 = 5 \times 13$ and $(-7)^2 + 8^2 = 113$, so the first could factor and the second does not.  
For $\mathfrak{p} = (a + bi)$ to check it is prime over $\mathbb{Z}[i]$, is it sufficient to check that $N(\mathfrak{p}) = a^2 + b^2$ is a prime over $\mathbb{Z}$?
It would be great to see the code in Sage or Pari/GP and that would be an acceptable answer.
 A: One could enumerate Gaussian integers within a range, say all $a + i b$
with $a$ and $b$ from $-30$ to $30$, and test whether they are prime and
satisfy the congruence (to be faster, enumerate only those that satisfy
the congruence, and test whether they are prime).
sage: G = GaussianIntegers()
sage: r = 30
sage: good = [] # prime and equal to 1 modulo (8 * ZZ[i])
sage: nope = [] # others
sage: for a in (-r .. r):
....:     for b in (-r .. r):
....:         g = G((a, b))
....:         if a % 8 ==1 and b % 8 == 0 and g.is_prime():
....:             good.append(g)
....:         else:
....:             nope.append(g)
....:

Then one could plot as follows:
sage: plot_good = point2d(good, color='red', size=10)
sage: plot_nope = point2d(nope, color='blue', size=1)
sage: (plot_good + plot_nope).show(aspect_ratio=1, figsize=7)
Launched png viewer for Graphics object consisting of 2 graphics primitives


or save the plot as a png file as follows:
sage: (plot_good + plot_nope).save('gaussian_primes.png', aspect_ratio=1, figsize=7)

Check our list for counterexamples to "prime in the Gaussian integers implies prime norm":
sage: [g for g in good if not g.norm().is_prime()]
[-23, -7]

The prime Gaussian integers $-7$ and $-23$ lie in
$1 + 8 \mathbb{Z}[i]$, but their norms are $7^2$ and $23^2$.
