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).
 A: 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. 
A: Denote by $\phi(D)$ the Jacobsthal sum
$$
\sum_t \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
\begin{align*}
\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).
\end{align*}
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
\begin{align*}
\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. 
\end{align*}
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).
$$
The sums in the brackets $\{ \}$ clearly are equal to $2(p-1)$; the last sum is $-2(p-1)$ because $\sum (\frac{st}p) = 0$, so
$$ \sum_{D}\phi(D)^2 = 2(p-1)^2 + 2(p-1) = 2(p-1)p. $$ 
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.)
A: (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).
A: 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.
