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

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.

-
Also: yes, the problem is not strictly comparable to (neither more or less general than) algorithms to solve the four-square problem. On the one hand, we require coprimality with $p$; on the other, we are not in characteristic zero. This does open the possibility of significantly different ways of obtaining solutions, while necessarily forcing new ones for "small" $n$. However, it is conceivable that an algorithm may be easily obtained by a modification of their algorithm. –  Niel de Beaudrap Mar 10 '10 at 21:56
Well: my approach to finding paths in the graphs $G_N$ is to experimentally try to decompose a residue into a sum or difference of squares, and decompose using the Chinese Remainder Theorem when I find factors of $N$. But using this reduction to find a path in $G_N$, seen as decomposition of a residue into quadratic units, I still require the same number of summands for each factor of $N$, which (as I remarked above) may be e.g. more than two. I picked 'four' for my question because I thought that solutions in $\mathbb N$ for the four-squares problem might be helpful. –  Niel de Beaudrap Mar 10 '10 at 22:41
Will: thanks for the references, they will be helpful. I still need to try and obtain an unconditional result, but I'm optimisitic now (esp. after re-reading Felipe's comments more carefully) that such an algorithm can be obtained. The problem is not that I have doubts as to whether the result will work 'empirically'; the problem is that I really want the algorithm, and not actually the paths for any particular case. Thus, I want a running time which is unconditional if such an algorithm can be obtained at all. –  Niel de Beaudrap Mar 13 '10 at 10:17

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

-

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?

-
Being a bit of a conservative sort, I would like an algorithm which works independently of unproven conjectures; although the idea that such a simple approach "ought to" work is a nice one. –  Niel de Beaudrap Mar 10 '10 at 21:59
Also: testing whether an integer is a quadratic residue currently doesn't have known efficient algorithms. (If it did, this problem would be substantially easier.) Not being familiar with the density-of-primes result which I suspect you might be using --- I am not myself a number theorist, and came upon this problem as if by accident --- does this amount to an attempt to solve the four-squares problem in the integers? (What, then, about e.g. the cases $n \in \{1,2,3\}$?) –  Niel de Beaudrap Mar 10 '10 at 22:07
Testing if a number $a$ is a quadratic residue mod p takes deterministically $O(\log p)$ steps by the stupid method of computing $a^{(p-1)/2}$ mod p (by square-and-multiply). Extracting the square root however can only be done in probabilistic polynomial time. My suggestion is not to solve the 4 squares problem for integers. You can even solve $n=x^2+y^2$ mod p for any nonzero n as long as p > 5. –  Felipe Voloch Mar 10 '10 at 22:18
Right. Quadratic residuacity is a difficult problem modulo composite N; the fact that the prime case is different slipped my mind. That just leaves the question of doing it unconditionally. –  Niel de Beaudrap Mar 10 '10 at 22:31
After reading your remarks more carefully (esp. with respect to root extraction), it seems clear to me that there is an efficient, unconditional, algorithm implicit in your responses. Thanks for reminding me that residuacity is easy modulo primes, and for the further information that root extraction is also easy with randomness modulo primes. –  Niel de Beaudrap Mar 13 '10 at 10:10

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.)

-