Let $p$ be an odd prime. An old theorem of Jacobi asserts that $p$ has exactly $8(p+1)$ representations as a sum of four squares of integers (solutions counted with order and sign). What is the most effective way to enumerate these solutions computationally? Can it be done in time $p^{1+\varepsilon}$, or even in time $p (\log{p})^A$?



Enumerating representations of an integer as a sum of squares
Form the set $S$ of all squares less than $p$. This has $O(\sqrt{p})$ elements, and writing them down takes $O(\sqrt{p} \log p)$ time. (You don't have to implement fast multiplication to do this; just compute the list of squares by successively adding odd numbers.)


Let $T$ be the set of all integers expressible as the sum of two elements of $S$. This has $O(p)$ elements, and takes $O(p \log p)$ steps to write down.


Sort $T$ and sort $p-T$. This is $O(p \log p)$ steps each. Find all duplicates between the lists $T$ and $p-T$; this takes $O(p)$ steps because they are already sorted. 


All in all, $O(p \log p)$ steps, the same size as the output.


