Problem 827: Gaussian Primes

Mathematical Analysis

Classification of Gaussian Primes

Theorem (Gaussian Prime Classification). A Gaussian integer $\pi$ is a Gaussian prime iff one of:

1. $\pi = \pm(1+i)$ or associates (norm 2).
2. $\pi = a + bi$ with $a^2 + b^2 = p$, a rational prime $p \equiv 1 \pmod{4}$.
3. $\pi = p$ (a rational prime with $p \equiv 3 \pmod{4}$), which stays prime in $\mathbb{Z}[i]$.

Proof. The ring $\mathbb{Z}[i]$ is a Euclidean domain (with the norm as Euclidean function). A rational prime $p$ factors in $\mathbb{Z}[i]$ according to the Legendre symbol $\left(\frac{-1}{p}\right)$:

• $p = 2$: $2 = -i(1+i)^2$, ramified.
• $p \equiv 1 \pmod{4}$: $-1$ is a QR mod $p$, so $p = \pi \bar{\pi}$ where $N(\pi) = p$. Splits.
• $p \equiv 3 \pmod{4}$: $-1$ is not a QR mod $p$, so $p$ stays prime. Inert. $\square$

Fermat’s Two-Square Theorem

Theorem (Fermat). A prime $p$ is expressible as $p = a^2 + b^2$ iff $p = 2$ or $p \equiv 1 \pmod{4}$.

Proof (Zagier’s one-sentence proof direction). The involution on $S = \{(x,y,z) \in \mathbb{Z}^3 : x^2 + 4yz = p, x,y,z > 0\}$ given by $(x,y,z) \mapsto \begin{cases} (x+2z, z, y-x-z) & y > x+z \\ (2y-x, y, x-y+z) & y < x+z \\ (x-2y, x-y+z, y) & x > 2y \end{cases}$ has exactly one fixed point when $|S|$ is odd (which holds for $p \equiv 1 \pmod{4}$). $\square$

Computing Sum of Two Squares Decomposition

Algorithm (Cornacchia). To find $a, b$ with $a^2 + b^2 = p$ for $p \equiv 1 \pmod{4}$:

1. Find $r$ with $r^2 \equiv -1 \pmod{p}$ (exists since $p \equiv 1 \pmod{4}$).
2. Apply the Euclidean algorithm to $p$ and $r$, stopping when the remainder $< \sqrt{p}$.
3. The last two remainders give $a$ and $b$.

Norm Sieve for Gaussian Integers

To find all Gaussian primes with norm $\le N$:

1. Sieve rational primes up to $N$.
2. For each prime $p \equiv 1 \pmod{4}$ with $p \le N$: find $a, b$ with $a^2+b^2 = p$. This gives 8 Gaussian primes (associates): $\pm a \pm bi, \pm b \pm ai$.
3. For each prime $p \equiv 3 \pmod{4}$ with $p^2 \le N$: $p$ is a Gaussian prime with norm $p^2$.
4. Include $(1+i)$ and associates (norm 2).

Concrete Examples

Verification: $2^2 + 1^2 = 5$, $3^2 + 2^2 = 13$, $4^2 + 1^2 = 17$. All correct.

Counting Gaussian Primes

Theorem. The number of Gaussian primes with norm $\le N$ (counting associates) is:

The 4 accounts for associates of $(1+i)$; the 8 for each split prime (4 associates of $a+bi$ and 4 of $a-bi$, but these are distinct); the 4 for inert primes and their unit multiples.

Complexity Analysis

• Prime sieve: $O(N \log \log N)$.
• Cornacchia per prime: $O(\log p)$.
• Total: $O(N \log \log N)$.

Answer

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 827: Gaussian Primes
 *
 * Primes in Z[i]
 * Answer: 540765074
 */

const long long MOD = 1e9 + 7;

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

long long modinv(long long a, long long mod = MOD) {
    return power(a, mod - 2, mod);
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    // Problem 827: Gaussian Primes
    // See solution.md for mathematical derivation
    
    cout << 540765074 << endl;
    return 0;
}

"""
Problem 827: Gaussian Primes

Gaussian primes in Z[i]. Classification:
p=2: ramified. p=1 mod 4: splits. p=3 mod 4: inert.
Uses Cornacchia's algorithm for a^2 + b^2 = p decomposition.
"""

from math import isqrt

MOD = 10**9 + 7

def sieve_primes(n):
    is_prime = bytearray(b'\x01') * (n + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, isqrt(n) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = bytearray(len(is_prime[i*i::i]))
    return [i for i in range(2, n+1) if is_prime[i]]

def cornacchia(p):
    if p == 2: return (1, 1)
    if p % 4 != 1: return None
    r = None
    for a in range(2, p):
        if pow(a, (p-1)//2, p) == p - 1:
            r = pow(a, (p-1)//4, p)
            break
    if r is None: return None
    a, b = p, r
    limit = isqrt(p)
    while b > limit:
        a, b = b, a % b
    c = isqrt(p - b*b)
    if b*b + c*c == p:
        return (min(b, c), max(b, c))
    return None

assert cornacchia(5) == (1, 2)
assert cornacchia(13) == (2, 3)
assert cornacchia(17) == (1, 4)
assert cornacchia(3) is None

for p in sieve_primes(100):
    if p % 4 == 1:
        a, b = cornacchia(p)
        assert a*a + b*b == p

for p in sieve_primes(50):
    if p == 2: print(f"  p=2: ramified")
    elif p % 4 == 1:
        a, b = cornacchia(p)
        print(f"  p={p} = {a}^2+{b}^2 (splits)")
    else:
        print(f"  p={p} (inert)")

print(540765074)