Lagrange's Theorem

Date: 02/27/2001 at 06:01:25
From: Andrew Chapman
Subject: Lagrange's theorem

Hi,

I see in your archives you show the proofs of Lagrange's theorem that 
every positive integer can be expressed as the sum of four squares - 
but is there an algorithm or some useful method for identifying 
*which* four squares? (I'm trying to write a computer program, you 
see...)

Thanks for your help!

Andrew

Date: 02/27/2001 at 16:07:54
From: Doctor Rob
Subject: Re: Lagrange's theorem

Thanks for writing to Ask Dr. Math, Andrew.

Indeed there is. In fact, there are at least three that I know.

Suppose you want to write N as a sum of four squares.

If N = 0 (mod 4), then write N/4 = a^2 + b^2 + c^2 + d^2, and
then N = (2*a)^2 + (2*b)^2 + (2*c)^2 + (2*d)^2.

If N = 2 (mod 4), then write N/2 = a^2 + b^2 + c^2 + d^2, and
then N = (a+b)^2 + (a-b)^2 + (c+d)^2 + (c-d)^2.

If N = 3 (mod 4), then let

   a = 2*[sqrt(N/4)] - 1,
   b = 2*[sqrt((N-a^2)/4)] - 1,
   M = N - a^2 - b^2.

You have reduced the problem to one of writing M = 1 (mod 4) as the
sum of two squares (see below). If this turns out to be impossible,
subtract 2 from b, recompute M, and try again, until you are able
to do this.

If N = 1 (mod 4), then let

   a = 2*[sqrt(N/4)],
   b = 2*[sqrt((N-a^2)/4)],
   M = N - a^2 - b^2.

Again you have reduced the problem to one of writing M = 1 (mod 4)
as the sum of two squares (see below).  If this turns out to be
impossible, subtract 2 from b, recompute M, and try again, until you
are able to do this.

There are a couple of ways to write M = 1 (mod 4) as the sum of two
squares. One is to factor M. If every prime factor of M that is
congruent to 3 (mod 4) appears to an even power, this is possible.
If any appears to an odd power, this is impossible.  Then you can
write M = u^2*P, where P is squarefree and has all its prime factors
congruent to 1 (mod 4). If you can write P as the sum of two squares, 
you are done. If you can write each prime factor of P as the sum of 
two squares, then you can use the identity

   (w^2 + x^2)*(y^2 + z^2) = (w*y+x*z)^2 + (w*z-x*y)^2

repeatedly to write P as the sum of two squares. That reduces the
problem to that of writing a prime number congruent to 1 (mod 4) as
the sum of two squares.

To solve

   p = a^2 + b^2, where p is prime, p = 1 (mod 4),

first solve

   x^2 = -1 (mod p).

Then start the Euclidean Algorithm with inputs p and x.  a and b
will be the first two remainders less than sqrt(p). That reduces the 
problem to solving x^2 = -1 (mod p).

You can do this by finding by trial and error a number y such that
the Legendre symbol (y/p) = -1.  Then x = y^((p-1)/2) (mod p) will
be such a solution. Since half the numbers y from 1 to p-1 will have 
this property, it is quick to find one.

A second method is to compose two lists, one of x^2 (mod N), and the 
other of -1-x^2 (mod N). Make the lists long enough so that there has 
to be a common element. That will be so if they cover the range 
0 < x <= N/2 + 1. Sort both lists, and find a common element. Those 
values of x[1] and x[2] will give you

   x[1]^2 + x[2]^2 + 1 = 0 (mod N),
   m*N = x[1]^2 + x[2]^2 + 1^2 + 0^2,

for some positive integer m. Thus m*N has been expressed as a sum of 
four squares. The trick is to replace m with smaller and smaller 
multiples until you get m = 1.

If m is even, then exactly one of x[1] and x[2] is even. Suppose
x[1] is even. Then define

   y[1] = x[1]/2,
   y[2] = x[1]/2,
   y[3] = (x[2]+1)/2,
   y[4] = (x[2]-1)/2.

It is easy to see that now

   (m/2)*N = y[1]^2 + y[2]^2 + y[3]^2 + y[4]^2,

and we have reduced m to a smaller number.

If m is odd, and m >= 3, define

   y[i] = x[i] (mod m), -(m-1)/2 <= y[i] <= (m-1)/2.

Then

   0 = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 (mod m),
     = y[1]^2 + y[2]^2 + y[3]^2 + y[4]^2 (mod m),
   m*m' = y[1]^2 + y[2]^2 + y[3]^2 + y[4]^2

for some integer m', and

   m*m' <= 4*(m-1)^2/4 = (m-1)^2,
   m' <= (m-1)^2/m < m - 1.

Also m' > 0.  Now

   (m*m')*(m*N) = (y[1]^2 + y[2]^2 + y[3]^2 + y[4]^2)*
                    (x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2)
                = z[1]^2 + z[2]^2 + z[3]^2 + z[4]^2,

where

   z[1] = (x[1]*y[1] + x[2]*y[2] + x[3]*y[3] + x[4]*y[4])/m,
   z[2] = (x[1]*y[2] - x[2]*y[1] + x[3]*y[4] - x[4]*y[3])/m,
   z[3] = (x[1]*y[3] + x[2]*y[4] - x[3]*y[1] - x[4]*y[2])/m,
   z[4] = (x[1]*y[4] + x[2]*y[3] - x[3]*y[2] - x[4]*y[1])/m.

It can be easily seen that the z[i]'s are integers.  Then

   m'*N = z[1]^2 + z[2]^2 + z[3]^2 + z[4]^2,

and m' < m, so we have reduced m.  Continue this until you reach
m = 1, and you're done.

There is another method that avoids factoring M, and also avoids
constructing and sorting the two long lists, if either is too 
difficult. If you need to know it, write again.

By the way, I believe the fact that every positive integer can be
written as the sum of four squares is called Bachet's Theorem,
since he first stated it explicitly in 1621. Lagrange supplied the 
first proof about a century and a half later.

- Doctor Rob, The Math Forum
  http://mathforum.org/dr.math/