Problem 971: Gaussian Integer Primes

Mathematical Foundation

Theorem 1 (Fermat’s Two-Square Theorem). An odd rational prime $p$ is expressible as a sum of two squares if and only if $p \equiv 1 \pmod{4}$.

Proof. ($\Leftarrow$) If $p \equiv 1 \pmod{4}$, then $-1$ is a quadratic residue modulo $p$ by Euler’s criterion: $(-1)^{(p-1)/2} \equiv 1 \pmod{p}$. Hence there exists $m$ with $m^2 \equiv -1 \pmod{p}$. Consider the lattice $\Lambda = \{(x,y) \in \mathbb{Z}^2 : y \equiv mx \pmod{p}\}$. This has index $p$ in $\mathbb{Z}^2$, so by Minkowski’s theorem applied to the disk of area $2\pi p > 4p = 4 \cdot \mathrm{vol}(\mathbb{Z}^2/\Lambda)$, there is a nonzero lattice point $(x,y) \in \Lambda$ with $x^2 + y^2 < 2p$. Since $y \equiv mx \pmod{p}$ implies $x^2 + y^2 \equiv x^2(1 + m^2) \equiv 0 \pmod{p}$, and $0 < x^2 + y^2 < 2p$, we get $x^2 + y^2 = p$.

($\Rightarrow$) If $p = a^2 + b^2$ with $p$ odd, then $a^2 + b^2 \equiv 0 \pmod{p}$. Since $\gcd(b,p)=1$, we get $(ab^{-1})^2 \equiv -1 \pmod{p}$, so $-1$ is a QR mod $p$, forcing $p \equiv 1 \pmod{4}$. $\square$

Theorem 2 (Classification of Gaussian Primes). A Gaussian integer $\pi = a + bi$ (up to units $\{\pm 1, \pm i\}$) is prime in $\mathbb{Z}[i]$ if and only if exactly one of the following holds:

1. $\pi = 1 + i$ (up to units and associates), corresponding to the ramified prime $2 = -i(1+i)^2$.
2. $a \neq 0$, $b \neq 0$, and $N(\pi) = a^2 + b^2$ is a rational prime (necessarily $\equiv 1 \pmod{4}$).
3. One of $a, b$ equals zero and the other is a rational prime $p \equiv 3 \pmod{4}$ (these are the inert primes).

Proof. The norm $N : \mathbb{Z}[i] \to \mathbb{Z}_{\ge 0}$ defined by $N(a+bi) = a^2 + b^2$ is multiplicative: $N(\alpha\beta) = N(\alpha)N(\beta)$. Since $\mathbb{Z}[i]$ is a Euclidean domain (with the norm as Euclidean function), it is a PID and hence a UFD.

Case (2): If $N(\pi) = a^2 + b^2 = p$ is a rational prime, and $\pi = \alpha\beta$ in $\mathbb{Z}[i]$, then $p = N(\alpha)N(\beta)$. Since $p$ is prime in $\mathbb{Z}$, one of $N(\alpha), N(\beta)$ equals 1, making that factor a unit. Hence $\pi$ is Gaussian prime.

Case (3): If $p \equiv 3 \pmod{4}$ is a rational prime and $p = \alpha\beta$ in $\mathbb{Z}[i]$, then $p^2 = N(p) = N(\alpha)N(\beta)$. If neither is a unit, $N(\alpha) = N(\beta) = p$, meaning $p$ is a sum of two squares, contradicting Theorem 1. So $p$ remains prime.

Completeness: Every Gaussian prime $\pi$ divides some rational prime $p$ (since $\pi \mid N(\pi) = \pi\bar{\pi}$ and $N(\pi)$ factors into rational primes). The splitting behavior $p = 2$: ramified; $p \equiv 1 \pmod{4}$: splits as $\pi\bar{\pi}$; $p \equiv 3 \pmod{4}$: inert --- exhausts all possibilities. $\square$

Lemma 1 (Counting in the First-Quadrant Square). For the region $\{a+bi : 0 \le a, b \le N\}$, the Gaussian primes consist of:

• Interior primes ($a > 0, b > 0$): pairs $(a,b)$ with $a^2 + b^2$ a rational prime.
• Axis primes on $b = 0$ ($a > 0$): rational primes $a \le N$ with $a \equiv 3 \pmod{4}$.
• Axis primes on $a = 0$ ($b > 0$): rational primes $b \le N$ with $b \equiv 3 \pmod{4}$.

Proof. This follows directly from Theorem 2, restricting the classification to the given region. Note that $1+i$ has $a = b = 1 > 0$ and $N(1+i) = 2$, which is prime, so it is counted among the interior primes. The point $(0,0)$ is not prime (it is zero). $\square$

Editorial

Count Gaussian primes a + bi in the first quadrant with 0 <= a, b <= 500. A Gaussian integer a + bi is prime if:. We sieve of Eratosthenes up to 2*N^2. Finally, interior Gaussian primes (a > 0, b > 0).

Pseudocode

Sieve of Eratosthenes up to 2*N^2
Interior Gaussian primes (a > 0, b > 0)
Axis primes on b = 0
Axis primes on a = 0

Complexity Analysis

• Time: $O(M \log\log M + N^2)$ where $M = 2N^2$. The sieve runs in $O(M \log\log M)$ and the double loop over the grid is $O(N^2)$. For $N = 500$, $M = 500{,}000$.
• Space: $O(M)$ for the prime sieve boolean array.

Answer

#include <bits/stdc++.h>
using namespace std;
int main() {
    const int N = 500, M = N*N*2;
    vector<bool> is_p(M+1, true); is_p[0]=is_p[1]=false;
    for(int i=2;(long long)i*i<=M;i++) if(is_p[i]) for(int j=i*i;j<=M;j+=i) is_p[j]=false;
    int cnt = 0;
    for(int a=1;a<=N;a++) for(int b=1;b<=N;b++) if(is_p[a*a+b*b]) cnt++;
    for(int a=1;a<=N;a++) if(is_p[a] && a%4==3) cnt += 2;
    cout << cnt << endl;
    return 0;
}

"""
Problem 971: Gaussian Integer Primes

Count Gaussian primes a + bi in the first quadrant with 0 <= a, b <= 500.

A Gaussian integer a + bi is prime if:
    - a > 0, b > 0 and a^2 + b^2 is a rational prime, OR
    - a > 0, b = 0 and a is a rational prime with a ≡ 3 (mod 4), OR
    - a = 0, b > 0 and b is a rational prime with b ≡ 3 (mod 4).

Key results:
    - Interior primes come from norms a^2 + b^2 being prime
    - Axis primes are primes p ≡ 3 (mod 4) placed on either axis
    - answer = total count in [0, 500] x [0, 500]

Methods:
    1. sieve_of_eratosthenes  — standard sieve up to needed bound
    2. count_interior_primes  — pairs (a,b) with a,b > 0 and a^2+b^2 prime
    3. count_axis_primes      — primes p ≡ 3 (mod 4) on the axes
    4. norm_distribution      — distribution of norms a^2+b^2 for Gaussian primes
"""

from math import isqrt
from collections import Counter

def sieve_of_eratosthenes(n):
    """Return bytearray where s[i] is 1 iff i is prime, for 0 <= i <= n."""
    s = bytearray(b'\x01') * (n + 1)
    s[0] = s[1] = 0
    for i in range(2, isqrt(n) + 1):
        if s[i]:
            s[i * i :: i] = bytearray(len(s[i * i :: i]))
    return s

def count_interior_primes(N, is_p):
    """Count pairs (a,b) with 1<=a,b<=N and a^2+b^2 prime. Also collect small ones."""
    count = 0
    small_a, small_b = [], []
    for a in range(1, N + 1):
        for b in range(1, N + 1):
            if is_p[a * a + b * b]:
                count += 1
                if a <= 80 and b <= 80:
                    small_a.append(a)
                    small_b.append(b)
    return count, small_a, small_b

def count_axis_primes(N, is_p):
    """Count Gaussian primes on both axes: (a,0) and (0,b) with a,b ≡ 3 mod 4 and prime."""
    count = 0
    axis_vals = []
    for a in range(1, N + 1):
        if is_p[a] and a % 4 == 3:
            count += 2  # (a, 0) and (0, a)
            axis_vals.append(a)
    return count, axis_vals

def norm_distribution(N, is_p, num_bins=50):
    """Histogram of norms a^2+b^2 for interior Gaussian primes."""
    norms = []
    for a in range(1, N + 1):
        for b in range(1, N + 1):
            n = a * a + b * b
            if is_p[n]:
                norms.append(n)
    return norms

#  Verification
# Small case: Gaussian primes with 0 <= a,b <= 5
_is_p = sieve_of_eratosthenes(50)
_int_count = 0
for a in range(1, 6):
    for b in range(1, 6):
        if _is_p[a * a + b * b]:
            _int_count += 1
_axis_count = 0
for a in range(1, 6):
    if _is_p[a] and a % 4 == 3:
        _axis_count += 2
_total_5 = _int_count + _axis_count
# Interior: (1,1)->2P (1,2)->5P (2,1)->5P (1,4)->17P (4,1)->17P (2,3)->13P (3,2)->13P
# (3,4)->25=5*5 no, (4,3)->25 no, (1,2)->5 yes, etc.
# Axis: 3 is prime and 3%4==3 -> +2
assert _axis_count == 2, f"Expected 2 axis primes for N=5, got {_axis_count}"

#  Main computation
N = 500
M = N * N + N * N  # max possible norm
is_p = sieve_of_eratosthenes(M)

interior_count, gp_a, gp_b = count_interior_primes(N, is_p)
axis_count, axis_vals = count_axis_primes(N, is_p)
answer = interior_count + axis_count
print(answer)

#  Visualization — 4-panel figure