A method for factoring integers into gaussian factors Is there a method for factoring a prime of the form $4k+1$ (residual 1 modulo 4) into Gaussian prime factors?
What I am looking for is a method/procedure to generate the factors from this table for the norm values which are, once again, primes of the form $4k+1$ (5, 13, 17, 29, ...)
There does not seem to be much on math stack exchange beyond this post, whose only comment is somewhat terse and I hope to see expanded upon.
Ex.
First few primes of the form $4k+1$ = $\{5, 13, 17, 29, ...\}$.
So for $5$, I want to know if there is a procedure for factoring it into $(2+i)(1+2i)$.
For $13$, how to factor it into $(3+2i)(2+3i)$, so on and so forth.
 A: Fermat's theorem on sums of two squares guarantees that every prime $ p \equiv 1 \pmod 4$ can be written as $p = a^2 + b ^2$. Moreover, this representation is also unique for positive $a,b$. This means that $p$ can be decomposed into Gaussian factors $p=(a+bi)(a-bi)$ uniquely (upto multiplication by units).
We know from quadratic reciprocity (or even Euler's criterion) that $-1$ is a square mod $p$ i.e. there exists an integer $t$ such that
$$ t^2 \equiv -1 \pmod p$$
The usual algorithm to find a square root mod $p$ is the Tonelli-Shanks algorithm. We only need a much simpler version in our case:
Take a random residue $a \pmod p$ and find $a^{\frac{p-1}{4}}$. Since $a^{\frac{p-1}{2}} \equiv \left( \frac{a}{p} \right) \pmod p$, we get a square root of $-1$ iff $a$ is a quadratic nonresidue modulo $p$. Since half of the residues mod $p$ are quadratic nonresidues, our average required number of tries is less than $2$.
Once we have found $t$, we see that $p \mid t^2+1=(t+i)(t-i)$. $p$ cannot completely divide one of $t \pm i$ (it has to divide all the coefficients), which means that each of the two factors of $p$ divide exactly one of $t \pm i$. WLOG assume $a+bi$ divides $t+i$ and since $\Bbb{Z}[i]$ is an Euclidean domain, we can carry out the Euclidean algorithm to find out $\operatorname{gcd}(t+i,p)$ to find $a+bi$. We easily generate $a-bi$ and the factorization is complete.
