Fermat's Christmas theorem on sums of two squares with Gaussian integers Gaussian integers are the set:
$$\mathbb{Z}[\imath] =\{a+b\imath : a,b \in\mathbb{Z} \}$$
With norm:
$$\mathrm{N}(a+b\imath)=a^{2}+b^{2}.$$
It satisfies $\mathrm{N}(\alpha\cdot \beta )=\mathrm{N}(\alpha)\cdot\mathrm{N}( \beta )$
The units of $\mathbb{Z}[\imath]$ are  precisely: $1,-1,\imath,-\imath$
Result: If $p=4k+1$ ( p : prime), there are $x<p$  such that $p\mid(x^{2}+1)$

$Theorem$: If $p$ is a prime with
$p=4k+1$, then $p$ is an sum of two
squares.
Proof.
If $p=4k+1$ there are $x<p$ such that $p\mid (x^{2}+1)$, then  $p\mid(x+\imath)(x-\imath)$.

Note that, if $p\mid (x+\imath)$, there are $(a+b\imath)$ such that:
$$(x+\imath)=p(a+b\imath)=pa+pb\imath.$$
This implies that $x=pa$, this is imposible as $x\lt p$.
Therefore $p\nmid(x+\imath)$, also $p\nmid(x-\imath)$.
Then $p$ is not an Gaussian prime, so $p = \alpha \beta$ with $\mathrm{N}(\alpha)\gt 1$ and $\mathrm{N}(\beta)\gt 1$.
If $ \alpha=a+b\imath $ and $\beta=c+d\imath$, we get:
$$\mathrm{N}(p)=\mathrm{N}(\alpha \beta)=\mathrm{N}\alpha\cdot\mathrm{N}\beta,$$
which implies:
$$p^{2}=(a^{2}+b^{2})(c^{2}+d^{2}).$$
Note that the last equation is an integer, so $(a^{2}+b^{2})|p^{2}$.
By this last equation $(a^{2}+b^{2}),(c^{2}+d^{2})\neq 1,$ so therefore:
$$p=a^{2}+b^{2}$$
**

It's beautiful!  Does anyone have other
proofs  ?

**
 A: 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.
A: 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.
A: There is a nice proof using Minkowski's theorem which also proves the four-square theorem.  It is Proof #2 here.
