One of their methods for $n=4k+2$ is as follows:

Randomly select an even number $a$ and an odd number $b$ such that $a^2+b^2 < n$. Then, we hope $p=n-a^2-b^2$ is a prime (you can show there's about a $1/(A \log n \log \log n)$ chance of $p$ being prime ); $p$ is of the form $4r+1$, so if $p$ is prime there's a solution to $c^2+d^2=p$.

To find that, we try to solve $m^2+1 \equiv 0 \pmod{p}$; I'm actually going to describe a slightly different method. Select $x$ at random from $1$ to $p-1$; then, $x^{2r}=\pm 1$ depending on whether $x$ is a quadratic residue so calculate $x^r$ by repeated squaring, and with a $1/2$ chance if $p$ is prime (a smaller one if $p$ is composite) you'll find a valid $m$.

Given such an $m$, $m+i$ is a Gaussian integer with norm divisible by $p$ but smaller than $p^2$; use the Euclidean algorithm on the Gaussian integers with $p$ and get $c+di$ with norm $p$, and $a^2+b^2+c^2+d^2=n$.

For an odd number $n$, solve for $a^2+b^2+c^2+d^2=2n$; note that by mod $4$ considerations exactly two of $a,b,c,d$ must be odd, and without loss of generality assume $a,b$ are odd and $c,d$ are even. Then,$(\frac{1}{2}(a+b))^2+(\frac{1}{2}(a-b))^2+(\frac{1}{2}(c+d))^2+(\frac{1}{2}(c-d))^2=n.$

For $n$ a multiple of $4$, solve for $n/4$ recursively, and multiply all values by $4$.

5more comments