Collapsing prime ideals of the form $\langle p, x+ i\rangle$ to a principal ideal We know $\mathbb{Z}[i]$ is a PID.
Also we know that if $p\equiv 1(4)$, $p$ splits in $\mathbb{Z}[i]$ and it factors into a product of prime ideals $\langle p, x+i \rangle \langle p, x-i \rangle$ where $x^2 \equiv -1(p)$. My question is how do we express these prime factors as principal ideals?
 A: This corresponds to the expression of $p$ as two squares: $p = a^2+b^2=(a+bi)(a-bi)$.
Since prime factorization is unique, we have $(p,x\pm i)=(a \pm bi)$.
Reversely, this gives an algorithm to compute $a$ and $b$ as the real and imaginary part of the gcd of $x+i$ and $p$.
A: Expressing these as principal ideals boil down to finding an element of norm $p$.
$x+iy$ is of norm $p$ if and only if $x^2+y^2 = p$.
You can solve this in $O(p)$ operations, either by going along the curve $x^2+y^2 = p$ step by step, or (if you already have determined that the ideal you want is $\langle p, k+i\rangle$), by replacing $y$ with $x/k \pmod p$ (because modulo that ideal, $i = -k = -x/y \pmod p$).
When you find a solution, the ideal it corresponds is the one with $k = x/y \pmod p$, and you obtain the other ideal by taking the conjugate.
If you are keeping a table of all the primes, another way to proceed is to simply start from the element $k+i$ and divide by the primes that are in excess.
So you look at the prime factorisation of its norm $k^2+1$, for each prime $q \neq p$ dividing it, you divide by an element generating the prime ideal $\langle q, (k \mod q)+i\rangle$.
Actually this is a good procedure to find the gaussian primes not necessarily in order by varying $k$ from $1$ to $\infty$. Each $k+i$ can have at most $1$ new prime, that you find by dividing by all the previously found primes.
