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

$\begingroup$ One value of p at a time, or for a range of p? It seems to make quite a difference. $\endgroup$ – Charles Matthews Jun 2 '10 at 16:52
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.

$\begingroup$ One caution occurs to me: can I actually sort that list in $O(p \log p)$ time? Sure, it has $p$ items, but they have length $\log p$, so I am not sure that I can treat comparisons as constant time. This is a very theoretical concern, though. Until $p$ overflows the word length of your computer, don't even think about it. $\endgroup$ – David E Speyer Jun 2 '10 at 17:08

$\begingroup$ Dear David, nicely done. Resembles the computation that finally settled Waring's problem for biquadrates. You might take a look at mathoverflow.net/questions/17711 by a quantum computation guy who evidently never actually programs anything. Felipe Voloch and I never figured out, well, what the guy's problem was. $\endgroup$ – Will Jagy Jun 2 '10 at 17:17

$\begingroup$ @David: No, if you sort $T$ the algorithm will be $O(p\log^2 p)$ time. But, there is no need to sort in the first place. Instead, insert $T$ into a hashtable, and go over $pT$ looking for each element in the hashtable. This is the desired minimum $O(p \log p)$ time. (When $p$ is greater than 64 bits this algorithm is quite impossible...) $\endgroup$ – Dror Speiser Jun 2 '10 at 22:25

$\begingroup$ I'm not completely convinced. If $x$ and $y$ are two random $N$ bit integers, then the expected time to determine which is larger is still $O(1)$: the probability that we need to consult $k$ bits in order to distinguish them is $2^{k}$. That heuristically makes me think that you might be able to sort in less than $O(p (\log p)^2)$ steps here. $\endgroup$ – David E Speyer Jun 2 '10 at 22:38

$\begingroup$ I am not sure whether or not a hash table is a better idea in practice. I'm certainly willing to believe that it is. But $T$ will have roughly $p/\sqrt{\log p}$ distinct elements, if I remember correctly, so I would expect the collision rate to be high. If you meant this as an "in theory" statement, isn't there usually a problem with theoretical worst case analyses in hashing being awful? But I am genuinely ignorant here. $\endgroup$ – David E Speyer Jun 2 '10 at 22:43
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.