Computing "Small" Gaussian Primes

simple.py

from sympy import isprime
from itertools import ifilter, count

def find_b(a):
    """Returns smallest b such that a^2 + b^2 is prime."""
    f = lambda b: isprime(a**2 + b**2)
    g = ifilter(f, count())
    b = g.next()
    return b

# for each a = 1,...,N compute the smallest b
# such that a+bi is a Gaussian prime
N = 10000
a = range(1, N)
b = map(find_b, a)