The theorem I am referring to is the following:
If $m$ and $n$ are integers, and $p$ is a prime number such that $4 \mid p – 3$, then the exponent of $p$ in the prime decomposition of $m^2 + n^2$ is even.
I can prove this using Fermat's Little Theorem (FLT), but I am looking for a more elementary proof.
This is because, in a book on (very) elementary number theory1, the reader is asked to prove this theorem as an exercise, and by this point the book has not covered anything remotely like FLT (or any modular arithmetic, for that matter).
In fact, by this point, this book has only covered the Euclidean algorithm, Bézout's identity, the uniqueness of prime decomposition, and some basic facts about Gaussian integers (including the uniqueness of the prime decomposition of Gaussian integers). The book has also stated, without proof, that an odd prime integer $p$ can be expressed as the sum of two squares if and only if $4 \mid p – 1$.
For what it's worth, at the point where I would otherwise apply FLT, I have reduced the problem to an equality of the form
$$ a^2b = c^2 + d^2$$
…where, $a$, $b$, $c$, and $d$ are integers, $\gcd(c, d) = 1$, and $b$ is square-free, and I must show that no odd prime factor $p$ of $b$ can satisfy $4 \mid p – 3$.
1 Weissman's An illustrated theory of numbers.
Best Answer
Here I'll use only some of the results that are listed in the body of the question as being known, with the exception that I make implicit use of the result (surely also known?) that irreducibles in $\mathbb{Z}[i]$ are prime.
First, I give two proofs that if $p \equiv 3 \pmod4$ then $p$ is a prime in $\mathbb{Z}[i]$.
This is for no better reason than that I only thought of the shorter proof after writing out the longer one! :(
Both proofs begin with an arbitrary odd prime $p$.
The first proof can of course be skipped without loss.
Stupidly long proof
We know that $p$ is a product of primes in $\mathbb{Z}[i]$.
If $\pi$ is a prime factor of $p$ in $\mathbb{Z}[i]$, then so is $\overline{\pi}$.
If $\overline{\pi} = \pi$, i.e., if $\pi$ is real, then $\pi$ is a factor of $p$ in $\mathbb{Z}$; but $\pi \ne 1$, therefore $\pi = p$.
Therefore, either $p$ is a prime in $\mathbb{Z}[i]$, or else none of the prime factors of $p$ in $\mathbb{Z}[i]$ are real.
If $\pi$ is a non-real factor of $p$ in $\mathbb{Z}[i]$, then $\pi$ and $\overline{\pi}$ cannot be associates, because then they would both be associates of $1 + i$, and $p$ could only be $2$, contradicting the postulate that $p$ is odd.
Therefore, either $p$ is a prime in $\mathbb{Z}[i]$, or else $p$ is a product of the form: \begin{align*} \pi_1\pi_2\cdots\pi_n\overline{\pi_1}\overline{\pi_2}\cdots\overline{\pi_n} & = \alpha\overline{\alpha} \\ & = a^2 + b^2 \\ & \equiv 0 \text{ or } 1 \pmod4, \end{align*} where $\alpha = \pi_1\pi_2\cdots\pi_n = a + ib$.
Therefore, either $p$ is a prime in $\mathbb{Z}[i]$, or else $p \equiv 1 \pmod4$.
That is, if $p \equiv 3 \pmod4$, then $p$ is a prime in $\mathbb{Z}[i]$. $\square$
(We have not needed to use the non-trivial result that a prime $\equiv 1 \pmod4$ is a sum of two squares.)
More sensible proof
As usual, denote the norm of a Gaussian integer $\alpha = a + ib$ by: $$ \operatorname{N}(\alpha) = a^2 + b^2 \equiv 0 \text{ or } 1 \pmod4. $$
If $p = \alpha\beta$ in $\mathbb{Z}[i]$, and neither $\alpha$ nor $\beta$ is a unit, then: $$ p^2 = \operatorname{N}(p) = \operatorname{N}(\alpha)\operatorname{N}(\beta), $$ where $\operatorname{N}(\alpha) \ne 1$, $\operatorname{N}(\beta) \ne 1$. Therefore, $\operatorname{N}(\alpha) = \operatorname{N}(\beta) = p$. Therefore, $p \equiv 1 \pmod4$.
So, if $p \equiv 3 \pmod4$, then either $\alpha$ nor $\beta$ must be a unit. That is, $p$ is a prime in $\mathbb{Z}[i]$. $\square$
(I use the word "prime" in this way because the word "irreducible" is not used in the question.)
Main result
From now on, we assume that $p \equiv 3 \pmod4$.
Given integers $m$ and $n$, we have: $$ m + in = p^r\beta, $$ where $r$ is a non-negative integer, and $\beta$ is a product of primes $\ne p$ in $\mathbb{Z}[i]$.
Therefore, $m - in = p^r\overline{\beta}$, where $\overline{\beta}$ is also a product of primes $\ne p$ in $\mathbb{Z}[i]$.
Therefore, the exponent of $p$ in the prime factorisation of $m^2 + n^2$ in $\mathbb{Z}[i]$ is $2r$.
But if $m^2 + n^2 = p^{2r}h$, where $h = \beta\overline{\beta}$, we cannot have $p | h$ in $\mathbb{Z}$, because then either $p | \beta$ or $p | \overline{\beta}$ in $\mathbb{Z}[i]$, contradicting the definition of $\beta$.
Therefore, the exponent of $p$ in the prime factorisation of $m^2 + n^2$ in $\mathbb{Z}$ is also $2r$. $\square$