Fermat's Christmas theorem on sums of two squares with Gaussian integers

This is best viewed from a slightly more general perspective as follows. In any $\rm UFD$, if $\rm\ a\ $ is not prime, i.e. $\rm\ a\:|\:bc\ $ but $\rm\ a\nmid b,\ a\nmid c\ $ then $\rm\ \gcd(a,b)\ $ is a proper factor of $\rm\:a\:$. Moreover this gcd can be computed when a $\rm UFD$ has a constructive Euclidean algorithm, as does $\rm \mathbb Z[i]\:$. Therefore this yields a constructive proof that nonprime nonunits are reducible in Euclidean domains (i.e. the nontrivial half of the equivalence of irreducible and prime elements in $\rm UFDs$).

Applying this above we deduce that $\rm\ gcd(p,x-i)\ = a + bi $ is a proper factor of $\rm\:p\:$, so it must have norm a proper factor of $\rm\ N(p) = p^2\ $, i.e. it must have norm $\rm\:p\:$. Therefore $\rm\ p = a^2 + b^2\ $ as desired. This leads to an elegant efficient few-line Euclidean algorithm to compute a representation of a prime $\rm\ p = 4k+1\ $ as a sum of two squares - see my post here for an implementation.

There is a proof by Roger Heath-Brown and another by John Ewell.

Heath-Brown's proof was later adapted into a "one-sentence" proof by Zagier. Zagier's proof is available at the wikipedia link given in Sivaram's answer.

There is a nice proof using Minkowski's theorem which also proves the four-square theorem. It is Proof #2 here.