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$?

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 $pT$. This is $O(p \log p)$ steps each. Find all duplicates between the lists $T$ and $pT$; this takes $O(p)$ steps because they are already sorted. All in all, $O(p \log p)$ steps, the same size as the output. 


Subtract a square you've not seen before and check for form 4^k(8m + 7)? If not, you have a sum of three squares. If yes, this square occurs in no decomposition as sum of four squares. After a number of trials only as large as the square root of p, you seem to have a list of possible summands and a simpler such question. 

