Can we find Gaussian primes $\pi = 1 + 8 \mathbb{Z}[i]$ with $N(\pi) < 10000$?

edit

Let $a+bi$ be a prime in $\Bbb Z[i]$. We can multiply it with $\pm1$ to have $a\ge 0$. Also, we can pass to the conjugate to have also the case $b\ge 0$. After these norming steps, $0\le a,b$. In case $a < b$, one can multiply with $i$ and renorm. So we are searching for primes of the form $a+ib$ where $0\le b\le a$.

Warming-up: It is simple to see if a number in $\Bbb Z[i]$ is prime or not. Dialog with the sage interpreter:

sage: K.<j> = QuadraticField( -1 )
sage: K(2018 + 2019*j).factor()
(-1) * (86*j - 139) * (-6*j - 5) * (-j - 2)
sage: K(2017).factor()
(-9*j - 44) * (9*j - 44)
sage: 2020.next_prime()
2027
sage: K(2027).factor()
2027
sage: K(2027).is_prime()
True

sage: K(1+8*j).factor()
(-3*j - 2) * (-j - 2)
sage: K(1+8*j).is_prime()
False

sage: K(-7+8*j).factor()
(j) * (7*j + 8)
sage: K(-7+8*j).is_prime()
True

To collect now all "normed" primes of the given shape with "real" and "imaginary" part up to a given bound we can try:

K.<j> = QuadraticField( -1 )
BOUND = 40
my_primes = []

for a in xrange( BOUND/8 ):
    for b in xrange( 1+8*a+1 ):
        p = 1 + 8*a + b*j
        if p.is_prime():
            my_primes.append(p)

for p in my_primes:
    print p

Above, there is a "decent bound", so i can print here the results:

j + 1
4*j + 9
2*j + 17
8*j + 17
10*j + 17
12*j + 17
4*j + 25
6*j + 25
12*j + 25
14*j + 25
16*j + 25
22*j + 25
24*j + 25
2*j + 33
8*j + 33
20*j + 33
28*j + 33
32*j + 33