How to find Gaussian primes $10^6 < |\mathfrak{p}| < 2 \times 10^6$ with Pari-GP?

Question

I am looking through the PARI/GP reference card.

Example, listing primes over $\mathbb{Z}$ is standard computer homework exercise [1]:

• 100129
• 126691
• 193013
• 284933
• 492113
• 769591

How to write a program to find Gaussian prime numbers $\mathfrak{p} = a + bi \in \mathbb{Z}[i]$ with $10^6 < |a|, |b| < 2 \times 10^6$ ?

---

Related:

• Prime elements in $\mathbb{Z}[\sqrt{2}]$

• Prime elements of $\mathbb{Z}[i\sqrt5]$.

• Finding Gaussian Primes along Lines in $\mathbb Z[i]$

Answer 1

I will compute below the numbers $\mathfrak p=a+ib$ constrained to the following "triangle" restrictions: $$ 10^6< b<a<2\cdot 10^6\ .$$ (From here build conjugates and associates to get also $\pm a\pm ib$ and $\pm b\pm ia$.)

In particular, $\mathfrak p\not\in\Bbb Z$, so $p:=\mathfrak p\bar{\mathfrak p}=a^2+b^2$ is a prime in $\Bbb Z$ which splits in two parts in the Gaussian integers. Here, $p$ is congruent to $1$ modulo four lying in the interval $J$ between $2\cdot 10^{12}$ and $8\cdot 10^{12}$.

One (slow) algorithm would take all these primes $p$ one by one, write them in the form $p=a^2+b^2$ and check the "triangle" condition. There are roughly some

? floor(intnum(x=2*10^12, 8*10^12, 1/log(x)))
%3 = 205708205954

primes in $J$, and half of them split in $\Bbb Z[i]$. A lot of primes for my taste. Below, there is an old fashioned pari/gp code implementing this slow algorithm. (It avoids the prime $1+i$, just in case somebody wants to relax $b<a$ to $b\le a$...)

The "real job" is done by the pari/gp factorization of Gaussian integers, a prime p in $\Bbb Z$ was slightly changed to I*p to implicitly tell pari/gp to factor over $\Bbb Z[i]$.

{f(N) = local(p, z, a, b, c, count);
count = 0;
forprime(p=2*N^2, 8*N^2,
    if( (p-3) % 4,
        z = factor(I*p)[1, 1];
        a = real(z); b = imag(z);
        if(a < b, c = a; a = b; b = c;);
        if( (N < b) && (a < 2*N), 
            z = a + I*b; count = count + 1;
            print(count, " --> ", z, " with norm ", p);)))}

Then f(100) delivers some $579$ such (normed) primes, and there is a quick end. The last ones are:

570 --> 190 + 189*I with norm 71821
571 --> 194 + 185*I with norm 71861
572 --> 195 + 184*I with norm 71881
573 --> 197 + 182*I with norm 71933
574 --> 195 + 188*I with norm 73369
575 --> 199 + 186*I with norm 74197
576 --> 196 + 191*I with norm 74897
577 --> 196 + 195*I with norm 76441
578 --> 199 + 194*I with norm 77237
579 --> 199 + 196*I with norm 78017
? 

Seen from this perspective, this algorithm seems to be "better" than the brute force algorithm considering all $(a,b)$ with $100<b<a<200$, and testing for the primality of $a+ib$. At the cost of having a good primes generator for $\Bbb Z$.

I have also started...

? f(10^6)
1 --> 1000006 + 1000001*I with norm 2000014000037
2 --> 1000008 + 1000003*I with norm 2000022000073
3 --> 1000010 + 1000003*I with norm 2000026000109
4 --> 1000008 + 1000007*I with norm 2000030000113
5 --> 1000013 + 1000002*I with norm 2000030000173
6 --> 1000013 + 1000008*I with norm 2000042000233
7 --> 1000019 + 1000004*I with norm 2000046000377
8 --> 1000019 + 1000014*I with norm 2000066000557
9 --> 1000018 + 1000017*I with norm 2000070000613
10 --> 1000021 + 1000014*I with norm 2000070000637
11 --> 1000026 + 1000009*I with norm 2000070000757
12 --> 1000028 + 1000007*I with norm 2000070000833
13 --> 1000021 + 1000016*I with norm 2000074000697
14 --> 1000027 + 1000010*I with norm 2000074000829

and the primes are slowly printed one by one, as they are fished and tested for relevant...