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 algorithm of 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 this graph. 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$. 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.