Gaussian primes

Gaussian integers

i =sqrt(-1)
3+4*i.conjugate() # conjugate see also Complex Numbers
3+4*i.imag() # imaginary part
3+4*i.real() # real part
3+4*i.norm() # norm

R.<x> = QQ[]; # Q[x]
K.<i> = NumberField(x^2 +1);
K; i^2

List Gaussian primes

Here , we call a Gaussian integer alpha a Gaussian prime if it is a prime element in Z[sqrt{-1}] and is in the range 0 < Re(alpha) and 0 < Im(alpha) < Re(alpha). A Gaussian integer a + bi is prime if and only if (i) a = b = 1, (prime above 2), (ii) b = 0, a is a rational prime = 1 mod 3 (unramified), or (ii) Norm(a+bi) = a^2 + b^2 is a rational prime = 3 mod 4 (comp. decomposed). 

# list Gaussian prime alpha with | alpha | < K
K = 100
i = sqrt(-1)
prime = [n for n in [2,..,K] if is_prime(n)]
GP = [p + 0*i for p in prime if p%4==3]
GP.append(1+i)
for a in [2,3,..,K]:
    for b in [1,2,..,a-1]:
        N = a^2 + b^2
        if is_prime(N) and  N< K^2:
            GP.append(a+b*i)

list_plot(GP).show(xmin=0,ymin=0,xmax=K,ymax=K)
According to the def, Gaussian primes correspondence to rational primes uniquely 


# list Gaussian prime alpha with prime below < K
t= cputime(subprocesses=True)
K = 10^6
i = sqrt(-1)
P = prime_range(3,K)
GP = [1+i]
for p in P:
   if p%4 == 3:
       GP.append(p+0*i)
   else:
       for a in [ceil(sqrt(p/2))..floor(sqrt(p))]:
           b = p - a^2
           if is_square(b):
               GP.append(a+sqrt(b)*i)
               break


#print(GP);
print len(GP)   # = prime_pi(K) (= # of primes less than K)
cputime(subprocesses=True)-t
K = 10^6: 25.710231 sec (MBA 11 2013) 

I = sqrt(-1)
GP=[1+I]
N=10^4
for p in prime_range(3,N):
    if p%4==3:
       c.append(p)
    else:
        for a in [ceil(sqrt(p/2))..ceil(sqrt(p))]:
           b = p - a^2
           if is_square(b):
               GP.append(a+sqrt(b)*I)
               break

twin = []
for pi in GP: 
    if pi + (1+ I) in GP: 
        twin.append(pi)
        print pi ,"&", pi + (1+ I)
here, we call pi and pi + (1+ i) are twin if they are Gaussian primes