Lagrange four-squares theorem: efficient algorithm with units modulo a prime?

Question

I'm looking at algorithms to construct short paths in a particular Cayley graph defined in terms of quadratic residues. This has led me to consider a variant on Lagrange's four-squares theorem.

The Four Squares Theorem is simply that for any $n \in \mathbb N$, there exist $w,x,y,z \in \mathbb N$ such that $$ n = w^2 + x^2 + y^2 + z^2 . $$ Furthermore, using algorithms presented by Rabin and Shallit (which seem to be state-of-the-art), such decompositions of $n$ can be found in $\mathrm{O}(\log^4 n)$ random time, or about $\mathrm{O}(\log^2 n)$ random time if you don't mind depending on the ERH or allowing a finite but unknown number of instances with less-well-bounded running time.

I am considering a Cayley graph $G_N$ defined on the integers modulo $N$, where two residues are adjacent if their difference is a "quadratic unit" (a multiplicative unit which is also quadratic residue) or the negation of one (so that the graph is undirected). Paths starting at zero in this graph correspond to decompositions of residues as sums of squares.

It can be shown that four squares do not always suffice; for instance, consider $N = 24$, where $G_N$ is the 24-cycle, corresponding to the fact that 1 is the only quadratic unit mod 24. However, finding decompositions of residues into "squares" can be helpful in finding paths in the graphs $G_N$. The only caveat is that only squares which are relatively prime to the modulus are useable.

So, the question: let $p$ be prime, and $n \in \mathbb Z_p ( := \mathbb Z / p \mathbb Z)$. Under what conditions can we efficiently discover multiplicative units $w,x,y,z \in \mathbb Z_p^\ast$ such that $n = w^2 + x^2 + y^2 + z^2$? Is there a simple modification of Rabin and Shallit's algorithms which is helpful?

Edit: In retrospect, I should emphasize that my question is about efficiently finding such a decomposition, and for $p > 3$. Obviously for $p = 3$, only $n = 1$ has a solution. Less obviously, one may show that the equation is always solvable for $n \in \mathbb Z_p^\ast$, for any $p > 3$ prime.

Answer 1

There is also an unconditional deterministic polynomial-time algorithm to find $x,y,z,w \in \mathbf{F}_p^\times$ such that $x^2+y^2+z^2+w^2=n$, given any $n \in \mathbf{Z}$ and any prime $p \ge 7$.

First, given $a,b,c \in \mathbf{F}_p^\times$, Theorem 1.10 of Christiaan van de Woestijne's thesis lets one find an $\mathbf{F}_p$-point $P$ on the smooth conic $C \colon ax^2+by^2=cz^2$ in $\mathbf{P}^2$ over $\mathbf{F}_p$. The usual trick of drawing lines through $P$ and taking the second point of intersection with $C$ lets one parametrize $C(\mathbf{F}_p)$. At most $6$ points of $C(\mathbf{F}_p)$ have one of $x,y,z$ equal to $0$, so by trying at most $7$ lines, one finds a point on the affine curve $ax^2+by^2=c$ with $x,y \in \mathbf{F}_p^\times$.

Now, to solve $x^2+y^2+z^2+w^2=n$, choose $c \in \mathbf{F}_p \setminus \{0,n\}$, and apply the previous sentence to find $x,y,z,w \in \mathbf{F}_p^\times$ satisfying $x^2+y^2=c$ and $z^2+w^2=n-c$.

Answer 2

How about trying $x,y,z=1,\ldots,[\log p]^2$ (or some such bound) and testing if $n-x^2-y^2-z^2$ is a square modulo $p$? That should be efficient. Proving that it works might require GRH. Did you want an algorithm with a proof?

Answer 3

As indicated by Felipe (primarily in his responses to my comments of his solution above), the problem is actually easy modulo a prime $p > 3$. Here I outline an explicit random poly-time solution, depending on ideas contributed by him.

First, the special case $p = 5$. We can only express 0 as a sum of an even number of quadratic units (which in this case are $\pm 1$), and can only express the quadratic units themselves using exactly one or at least three quadratic units. We set this case aside and assume $p > 5$.

Second, for any $p \equiv 3 \pmod{4}$, the quadratic residues do not include $-1$; therefore $0$ cannot be formed as a quadratic residue. It suffices however to represent any $n \ne 0$ as a sum of two quadratic units, in which case we may easily reduce the problem to expressing $n \in \mathbb Z \setminus$ { 0 } as a sum of two quadratic units.

A classical result is that modulo $p$, and for $p > 5$, the number of ordered pairs $(q,q+1)$ of consecutive quadratic residues is asymptotically a large constant fraction of $p$ (specifically 1/4); and similarly for the number of pairs $(q,q+1)$ for which $q$ is a quadratic unit and $q+1$ a "non-quadratic" unit. We may then express $$ n = nq(q+1)^{-1} + n(q+1)^{-1} $$ where we choose $q+1$ to have the same "residuacity" as $n$; this is then a sum of two quadratic units. Because the suitable pairs $(q,q+1)$ occur a large constant fraction of the time, we may easily generate such a pair whether $n$ is a quadratic unit or not; and by efficiently finding roots for $nq(q+1)^{-1}$ and $n(q+1)^{-1}$, we may then find unconditionally find solutions in random polynomial time.

(I post this answer in order to provide an explicit record, and to emphasize that it can be done unconditionally; however I'm upvoting his answer as the one which contributed the useful ideas.)