An Explicit Solution to $a^2+b^2=p$ $\def \op {\operatorname*} \def \C#1#2 {\begin{pmatrix} #1\\#2 \end{pmatrix}}$
I got the following theory from the internet and I seek a proof:
Let $p=4k+1 (k \in \Bbb{Z}^+)$ be a given prime. Assume that $$2^{-1} \C{2k}{k} \equiv a \pmod{p},\quad -2k \le a \le 2k$$ then $\sqrt{p-a^2}$ is an integer. ($2^{-1}$ is seen as $\frac{p+1}{2}$, the inverse of $2$ in $\Bbb{Z}_p$)
I verified the statement for all $p<50000$ by programming. The question is that although $a,b$ seem to have equal status, this calculation can only find out one of the two numbers. Also, the calculation does not ensure that $a^2 \lt p$ but all the results meet this condition, so I suspect the combinatorial number to be an approximation of some transcendental function. I also hope to get some further studying if anyone learns about what is behind the method.
Bonus: A similar result is given when $p=6k+1 (k \in \Bbb{Z}^+)$ is a prime. Let $$2^{-1} \C{3k}{k} \equiv a \pmod{p},\quad -3k \le a \le 3k$$ then $\sqrt{\frac{p-a^2}{3}}$ is an integer. I think they should derive from the same part somewhere in number theory.
 A: As noted in comments, this can be shown using Jacobi sums. Let $\newcommand{\ffield}{\mathbb{F}_p}\ffield:=\mathbb{Z}/p\mathbb{Z}$; we compute
\begin{align*}
A_p&:=\#\{(x,y)\in\ffield^2:x^4+y^2=1\}&(p&=4k+1)
\\
B_p&:=\#\{(x,y)\in\ffield^2:x^3+y^2=1\}&(p&=6k+1)
\end{align*}
in two different ways. First, for $a\in\ffield$ we have $\#\{y\in\ffield:y^2=a\}=1+\chi_2(a)$, where $$\chi_2(a):=\Big(\frac{a}{p}\Big)\equiv a^{(p-1)/2}\pmod{p}$$ is the Legendre symbol. So, $A_p=\sum_{x\in\ffield}\big(1+\chi_2(1-x^4)\big)=p+\sum_{x\in\ffield}\chi_2(1-x^4)$, but $$\sum_{x\in\ffield}(1-x^4)^{2k}=\sum_{j=0}^{2k}\binom{2k}{j}(-1)^j\sum_{x\in\ffield}x^{4j}=(-1)^{k+1}\binom{2k}{k}-1$$ because $\sum_{x\in\ffield}x^n=0$ if $n$ is not a nonzero multiple of $p-1[=4k]$ (otherwise, the sum is $-1$; we assume $0^0=1$ everywhere), so that only the terms with $j=k$ and $j=2k$ of the outer sum "survive". Thus, $A_p\equiv(-1)^{k+1}\binom{2k}{k}-1\pmod{p}$. Similarly, for $B_p$ we get a sum with the only nonzero term at $j=2k$, ending up in $B_p\equiv-\binom{3k}{2k}=-\binom{3k}{k}\pmod{p}$.
Second, the Jacobi sums way. Here we recall that, for $a\in\ffield$ and $d\mid p-1$, we have $$\#\{x\in\ffield:x^d=a\}=\sum_{\chi^d=1}\chi(a),$$ with $\chi$ running over characters modulo $p$ whose order divides $d$ (this generalizes the case $d=2$ used above). So, for $p=4k+1$ and $\chi$ a character modulo $p$ of order $4$ (thus $\chi^2=\chi_2$), $$A_p=\sum_{\substack{a,b\in\ffield\\a+b=1}}\#\{x:x^4=a\}\cdot\#\{y:y^2=b\}=\sum_{\substack{0\leqslant r\leqslant 3\\0\leqslant s\leqslant 1}}J(\chi^r,\chi_2^s).$$
Computing nonzero terms using the properties of Jacobi sums, we have $$A_p=p+J(\chi,\chi_2)+J(\chi_2,\chi_2)+J(\chi^3,\chi_2)=p-1+2\Re J(\chi,\chi_2)$$ since $J(\chi_2,\chi_2)=-\chi_2(-1)=-1$ and $\chi^3=\bar\chi\implies J(\chi^3,\chi_2)=\overline{J(\chi,\chi_2)}$.
But $J(\chi,\chi_2)\in\mathbb{Z}[i]$ (because $\chi$ takes values from there) and $|J(\chi,\chi_2)|=p^{1/2}$ (known one, follows from the analogous property of the Gauss sum). That is, there are integers $a_p,b_p$ such that $J(\chi,\chi_2)=a_p+b_p i$ and $a_p^2+b_p^2=p$; the above gives $A_p=p-1+2a_p$. Combined with the "binomial" expression obtained earlier, and noting that $|a_p|=|a|$ for the "$a$" in the question, we're done.
Similarly, for $p=6k+1$ and $\chi$ a character modulo $p$ of order $3$, we get $$B_p=p+J(\chi,\chi_2)+J(\chi^2,\chi_2)=p+2\Re J(\chi,\chi_2),$$ and this time $J(\chi,\chi_2)\in\mathbb{Z}[\omega]$ with $\omega=e^{2\pi i/3}$. So, for some integers $a_p,b_p$, $$J(\chi,\chi_2)=a_p+b_p\omega,\qquad B_p=p+2a_p-b_p,\\4p=4|J(\chi,\chi_2)|^2=4(a_p^2-a_p b_p+b_p^2)=(2a_p-b_p)^2+3b_p^2,$$ and we obtain the Bonus claim (fixme: exclude the possibility of half-an-integer...).
