In a post Efficiently finding two squares which sum to a prime
I read
"In 1825 Gauss gave the following construction for writing a prime congruent to $1 \pmod{4}$ as a sum of two squares: Let $p=4k+1$ be a prime number. Determine $x$ (this is uniquely possible...) so that
$$ x = \frac{(2k)!}{2(k!)^2} \pmod{p}, \quad |x| < \frac{p}{2}$$
Now determine $y$ so that
$$ y = x \cdot (2k)! \pmod{p}, \quad |y| < \frac{p}{2}$$
Gauss showed that $x^2+y^2=p$."
I checked the Stark's book and I did not see a direct reference to the Gauss original paper or book or other reference to explain the method how to come to this and find the formula. If you know the reference or you can explain its procedure, I will be grateful.
We have $$ (2k)!^2 \equiv (-1)^{2k}(4k)! = (p-1)! \equiv -1 \pmod p $$ by writing $-1 \equiv 4k \pmod p$, $-2\equiv 4k-1 \pmod p$, etc. so that $(2k)! \equiv (-1)^{2k} (4k)(4k-1)\cdots(2k+1) \pmod p$. Then $$ x^2 \equiv \frac{(2k)!^2}{4(k!)^4} \equiv \frac{-1}{4(k!)^4} \pmod p. $$ Similarly, $y^2 \equiv \frac{1}{4(k!)^4} \pmod p$. So $x^2 + y^2 \equiv 0 \pmod p$. Now I am guessing that Gauss brilliantly bounded $x^2+y^2$ using his choice of $x$ and $y$ but the trivial bound is $$ x^2 + y^2 = |x|^2 + |y|^2 < \frac{p^2}4 + \frac{p^2}4 = \frac{p^2}2. $$ So I am guessing Gauss did something smart ; I checked the book referenced in your link and it doesn't redirect to a proof, the author only says that Gauss proved that.
Hope that helps,