Associated Topics || Dr. Math Home || Search Dr. Math

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

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/
```
Associated Topics:
College Number Theory
High School Number Theory

Search the Dr. Math Library:

 Find items containing (put spaces between keywords):   Click only once for faster results: [ Choose "whole words" when searching for a word like age.] all keywords, in any order at least one, that exact phrase parts of words whole words

Submit your own question to Dr. Math
Math Forum Home || Math Library || Quick Reference || Math Forum Search