Explicit formula for Fermat's 4k+1 theorem

Let $p$ be a prime number of the form $4k+1$. Fermat’s theorem asserts that $p$ is a sum of two squares, $p=x^2+y^2$.

There are different proofs of this statement (descent, Gaussian integers,…). And recently I’ve learned there is the following explicit formula (due to Gauss):
$x=\frac12\binom{2k}k\pmod p$,
$y=(2k)!x\pmod p$
($|x|,|y|<p/2$). But how to prove it?

Remark. In another thread Matt E also mentions a formula
$$x=\frac12\sum\limits_{t\in\mathbb F_p}\left(\frac{t^3-t}{p}\right).$$
Since $\left(\dfrac{t^3-t}p\right)=\left(\dfrac{t-t^{-1}}p\right)=(t-t^{-1})^{2k}\mod p$ (and $\sum t^i=0$ when $0<i<p-1$), this is, actually, the same formula (up to a sign).

Solutions Collecting From Web of "Explicit formula for Fermat's 4k+1 theorem"

Here is a high level proof. I assume it can be done in a more elementary way. Chapter 3 of Silverman’s Arithmetic of Elliptic Curves is a good reference for the ideas I am using.

Let $E$ be the elliptic curve $y^2 = x^3+x$. By a theorem of Weyl, the number of points on $E$ over $\mathbb{F}_p$ is $p- \alpha- \overline{\alpha}$ where $\alpha$ is an algebraic integer satisfying $\alpha \overline{\alpha} =p$, and the bar is complex conjugation. (If you count the point at $\infty$, then the formula should be $p – \alpha – \overline{\alpha} +1$.)

Let $p \equiv 1 \mod 4$. We will establish two key claims: Claim 1: $\alpha$ is of the form $a+bi$, for integers $a$ and $b$, and Claim 2: $-2a \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$. So $a^2+b^2 = p$ and $a \equiv -\frac{1}{2} \binom{(p-1)/2}{(p-1)/4}$, as desired.

Proof sketch of Claim 1: Let $R$ be the endomorphism ring of $E$ over $\mathbb{F}_p$. Let $j$ be a square root of $-1$ in $\mathbb{F}_p$. Two of the elements of $R$ are $F: (x,y) \mapsto (x^p, y^p)$ and $J: (x,y) \mapsto (-x,jy)$.

Note that $F$ and $J$ commute; this uses that $j^p = j$, which is true because $p \equiv 1 \mod 4$. So $F$ and $J$ generate a commutative subring of $R$. If you look at the list of possible endomorphism rings of elliptic curves, you’ll see that such a subring must be of rank $\leq 2$, and $J$ already generates a subring of rank $2$. (See section 3.3 in Silverman.) So $F$ is integral over the subring generated by $J$. That ring is $\mathbb{Z}[J]/\langle J^2=-1 \rangle$, which is integrally closed. So $F$ is in that ring, meaning $F = a+bJ$ for some integers $a$ and $b$.

If you understand the connection between Frobenius actions and points of $E$ over $\mathbb{F}_p$, this shows that $\alpha = a+bi$.

Proof sketch of Claim 2: The number of points on $E$ over $\mathbb{F}_p$ is congruent modulo $p$ to the coefficient of $x^{p-1}$ in $(x^3+x)^{(p-1)/2}$ (section 3.4 in Silverman). This coefficient is $\binom{(p-1)/2}{(p-1)/4}$. So
$$- \alpha – \overline{\alpha} \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$$
or
$$-2a \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$$
as desired.

Remark:
This is very related to the formula Matt E mentions. For $u \in \mathbb{F}_p$, the number of square roots of $u$ in $\mathbb{F}_p$ is $1+\left( \frac{u}{p} \right)$. So the number of points on $E$ is
$$p+\sum_{x \in \mathbb{F}_p} \left( \frac{x^3+x}{p} \right).$$

This is essentially Matt’s sum; if you want, you could use the elliptic curve $y^2 = x^3-x$ in order to make things exactly match, although that would introduce some signs in other places. So your remark gives another (morally, the same) proof of Claim 2.

There is a proof on page 192 of Franz Lemmermeyer’s book, Reciprocity Laws from Euler to Eisenstein. I found it by typing “sum of two squares” and “binomial coefficient” into Google.

There is also a proof in Allan Adler’s paper, Eisenstein and the Jacobian varieties of Fermat curves, which paper is freely available on the web.

Denote by $\phi(D)$ the Jacobsthal sum
$$\sum_D\left(\frac tp\right)\left(\frac{t^2+D}p\right).$$
Note that $\phi(a^2D)=(a/p)\phi(D)$, so $|\phi(D)|$ depends only on $\bigl(\frac Dp\bigr)$; let $2x$ and $2y$ be the absolute values of the Jacobsthal sum for quadratic residue and quadratic non-residue respectively.

Theorem (Jacobsthal, 1907). If $p=4k+1$, $x^2+y^2=p$.

(As explained in the OP, this statement implies the formula of Gauss.)

Proof. Obviously,
$$\sum_{D\in\mathbb F_p}\phi(D)^2=2(p-1)(x^2+y^2),$$
so we need to compute the sum
$$\sum_{D}\phi(D)^2= \sum_{D,s,t}\left(\frac{st}p\right)\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right)= \sum_{s,t}\left(\frac{st}p\right)\sum_D\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right).$$

Lemma.
$$\sum_D\left(\frac Dp\right)\left(\frac{D+c}p\right)=\begin{cases} \phantom p-1,&c\neq 0;\\ p-1,&c=0. \end{cases}$$

To prove the lemma for $c\neq0$ observe that in this case
$$\sum_D\left(\frac{D^2+cD}p\right)= \sum_{D\neq0}\left(\frac{1+cD^{-1}}p\right)= \sum_{D’\neq1}\left(\frac{D’}p\right)=-\left(\frac1p\right)=-1.$$

Now let us apply the lemma to our sum:
$$\sum_D\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right)= \sum_{D’}\left(\frac{D’}p\right)\left(\frac{D’+t^2-s^2}p\right)= \begin{cases} \phantom p-1,&s\neq\pm t;\\ p-1,&s=\pm t; \end{cases}$$
so
$$\sum_{D}\phi(D)^2=(p-1)\left\{\sum_{s=t}\left(\frac{t^2}p\right)+\sum_{s=-t}\left(\frac{-t^2}p\right)\right\}+(-1)\sum_{s\neq\pm t}\left(\frac{st}p\right).$$
For $p=4k+1$ this equals $2p(p-1)$ (the last sum being zero). Hence indeed $x^2+y^2=p$.

(To my surprise, although the proof is quite short and elementary, I haven’t been able to find a concise modern reference. Anyway, Jacobsthal’s paper «Über die Darstellung der Primzahlen der Form $4n+1$ als Summe zweier Quadrate» is quite accessible.)

(Below is a kind of a low-tech version of David Speyer’s answer — [essentially] proving Weil for a cubic instead of using it. On the other hand, it’s a kind of a more conceptual version of the computation with Jacobsthal sums.)

I. Preliminaries

For any two (nontrivial multiplicative) characters mod $p$ define Jacobi sum
$$J(\chi,\lambda):=\sum_{a+b=1}\chi(a)\lambda(b).$$
Recall that (if $\chi\lambda$ is a non-trivial character)
$$J(\chi,\lambda)=\frac{g(\chi)g(\lambda)}{g(\chi\lambda)},$$
where
$$g(\chi)=\sum\chi(t)\zeta^t$$
is the usual Gauss sum (cf. beta-function and gamma-function, btw).

Corollary. Since $|g(\chi)|=\sqrt p$, the absolute value of a Jacobi sum is $\sqrt p$ (given that $\chi$, $\lambda$ and $\chi\lambda$ are non-trivial).

II. Proof

Take $\chi$ to be a non-trivial character of order $4$ mod $p=4k+1$; $\chi^2=:\rho$ is then the quadratic character. Since $\chi$ and $\rho$ take values in $\{\pm1,\pm i\}$, the sum $J(\chi,\rho)$ is a Gaussian integer, and by Corollary it has norm $p$.

This already gives a proof of Fermat’s $4k+1$ theorem. And explicitly (using that $J$ is real only if $a$ is a quadratic residue)
$$\operatorname{Re}J(\chi,\rho)=\frac12\sum_t\left(\frac tp\right)\left(\frac{1-t^2}p\right)=\frac12\sum_t\left(\frac{t-t^3}p\right).$$
Which is, of course, just an example of counting of points on an elliptic curve ($y^2=x^3-x$ or $y^2=1-x^4$) — well, counting points on hypersurfaces is the point of Jacobi sums, after all.

Reference: Ireland, Rosen. A Classical Introduction to Modern Number Theory (ch. 8).