
Some of the four squares can be zero, so it's fair to say that the sum of 1, 2, or 3 squares can be expressed as the sum of 4 squares, of which some are zero.
The Quaternion Identity (see below) gives us a way to express the product of two numbers (each of which can be expressed as the sum of 4 squares) as the sum of 4 squares.
Since every number can be expressed as the product of primes, it's enough to show that every prime can be expressed as the sum of 4 squares.
This can be illustrated with the following practical method of expressing any number, n, as the sum of four squares:
First, factor out the largest square factor, s^{2}, leaving only distinct prime factors, p_{1}p_{2}p_{3}...p_{m}. (n=s^{2}p_{1}p_{2}p_{3}...p_{m})
Then use repeated applications of the Quaternion Identity. This is a cute trick in which two numbers, each of which are expressed as the sum of four squares, can be multiplied together and expressed as the sum of four squares. The hard part, which is the bulk of this web page, is finding four squares that add up to a given prime number. Once that's done for all the prime factors of n, then it's just a matter of multiplying them together, like this:
Let A_{0}=1, and B_{0}=C_{0}=D_{0}=0.
For each k, where k varies from 1 to m, do the following steps, which are repeated applications of the Quaternion Identity:
find four squares a_{k}^{2}, b_{k}^{2}, c_{k}^{2}, d_{k}^{2} whose sum is p_{k}.
Let A_{k}=(a_{k}A_{k1} + b_{k}B_{k1} + c_{k}C_{k1} + d_{k}D_{k1}),
Let B_{k}=(a_{k}B_{k1}  b_{k}A_{k1} + c_{k}D_{k1}  d_{k}C_{k1}),
Let C_{k}=(a_{k}C_{k1}  c_{k}A_{k1} + d_{k}B_{k1}  b_{k}D_{k1}),
Let D_{k}=(a_{k}D_{k1}  d_{k}A_{k1} + b_{k}C_{k1}  c_{k}B_{k1})Finally, the original number, n, can be expressed as the sum of four squares this way:
n = (sA_{m})^{2} + (sB_{m})^{2} + (sC_{m})^{2} + (sD_{m})^{2}
Now we know that any composite number can be expressed as the product of primes, so all we need to show is that any prime number can be expressed as the sum of four squares, and we have a method of expressing their products as the sum of four squares, too. Primes of the form 4k+1 can be expressed as the sum of two squares, as explained immediately below. Primes of the form 4k+3 can be expressed as the sum of three or four squares, as explained below that. That leaves only the smallest prime, 2, which can be expressed as 1^{2}+1^{2}.
Since any prime can be expressed as the sum of four squares, and any product of two sums of four squares can be expressed as the sum of four squares, our proof is not only complete, but constructive as well.
Use the procedure, below, to find a number, a, such that a^{2} = 1 (mod p).
Then, once you have found a square root of 1, mod p, then the following algorithm will give you a pair of squares whose sum is p.
a^{2} = 1 (mod p),
so a^{2}+1^{2}=kp, and since a < p, a^{2}+1^{2}<p^{2}, so k < p.The following is an iterative procedure which, given a^{2}+b^{2}=kp, finds a new pair of numbers A^{2}+B^{2}=k'p, where k' <= k/2. By repeating this procedure, we find values of k that are smaller and smaller, eventually reaching 1. In this way, we find a pair of numbers A and B such that A^{2}+B^{2}=p.
We start with a^{2}+b^{2}=kp, and call this equation (1).
Let A be the absolutely least residue of a mod k 
i.e. the smaller in absolute value of the residue a, mod k, and ak
let B be the absolutely least residue of b mod k.
Then A and B are in the interval (k/2,k/2], so
A^{2}+B^{2} <= (k/2)^{2}+(k/2)^{2} = k^{2}/2.
Also, since A=a (mod k) and B=b (mod k), A^{2}+B^{2}=a^{2}+b^{2}=0 (mod k), or in other words,
A^{2}+B^{2}=mk, with m<k, and call this equation (2).Multiplying the lefthand and righthand sides of equations (1) and (2),
(a^{2}+b^{2})(A^{2}+B^{2}) = mk^{2}p = (aA+bB)^{2}+(aBbA)^{2}, from the Complex Product Identity
aA+bB = a^{2}+b^{2} = 0 (mod k), and aBbA = abba = 0 (mod k), so we can divide through by k^{2}, giving us
((aA+bB)/k)^{2}+((aBbA)/k)^{2} = mpThis algorithm completes in log_{2}p time, resulting in a pair of squares adding up to p.
Example: find two numbers whose squares add up to 53
Start looking for an x such that x^{26}=1 (mod 53).
2^{26}=1 (mod 53), so let a=2^{13}=23 is a square root of 1, and let b=1...
a^{2} + b^{2} = kp. Filling in the numbers,
(23)^{2} + 1^{2} = 53k, so k=10
The absolutely least residue of 23 and 1 mod 10 are A=3 and B=1.
My new "a" and "b" are (aA+bB)/k and (aBbA)/k, which are
(23*(3)+1*1)/10=7 and (23*11*(3))/10=2.Starting over, with a=7 and b=2, we have a^{2} + b^{2} = kp:
7^{2} + (2)^{2} = 53k, so k=1. We've found our two numbers whose squares add up to 53.
We begin by finding a solution to y^{2}=1x^{2} (mod p). (See the proof that there's at least one solution, below.)
Empirically, we just start examining possible values for the righthandside: 12^{2}, 13^{2}, etc. until we find one that's a quadratic residue (i.e. a perfect square), mod p. In practice it seldom takes many tries, because about half of all these numbers are quadratic residues. And luckily, for primes of the form 4k+3, there's an excellent test to see if it's a quadratic residue.
From Fermat's Little Theorem, we know that
c^{p1} = 1 (mod p),
so c^{p+1} = c^{2} = (c)^{2} (mod p),
so c^{(p+1)/4} is a square root of c (mod p), or c^{(p+1)/4} is a square root of c (mod p).In practice, to see if a number, c, is a quadratic residue, mod p, we set b=c^{(p+1)/4}, and then check whether b^{2}=c. If it is, then we've found a solution to y^{2}=1x^{2}.
The test works by checking successively
Is (11^{2})^{(p+1)/2} = 11^{2} ?
Is (12^{2})^{(p+1)/2} = 12^{2} ?
Is (13^{2})^{(p+1)/2} = 13^{2} ?
...Until we find a solution such that (1x^{2})^{(p+1)/2} = 1x^{2}.
Then, since y^{2} = 1x^{2} = (1x^{2})^{(p+1)/2}, we know that y = (1x^{2})^{(p+1)/4}.We start by letting a be the absolutely least residue of y, b be the absolutely least residue of x, c=1, and d=0.
a < (p1)/2 and b < (p1)/2, so a^{2}+b^{2}+1 < p^{2}/2, so a^{2}+b^{2}+c^{2}+d^{2} = kp, with k<p.
Then we follow an iterative procedure similar to that for primes of the form 4k+1, except using the Quaternion Identity (see below) instead of the Complex Product Identity. In this procedure, given a^{2}+b^{2}+c^{2}+d^{2}=kp, we find a new set of four numbers A^{2}+B^{2}+C^{2}+D^{2}=k'p, where k' <= k/2.
If k is odd...
As before, we start with a^{2}+b^{2}+c^{2}+d^{2} = kp, and call this equation (1)
Let A, B, C, and D be the absolutely least residues mod k of a, b, c, and d.
Then each of A, B, C, and D is <= (k1)/2, so
A^{2}+B^{2}+C^{2}+D^{2} <= (k1)^{2}
A^{2}+B^{2}+C^{2}+D^{2} = a^{2}+b^{2}+c^{2}+d^{2} = 0 (mod k), so
A^{2}+B^{2}+C^{2}+D^{2} = mk, with m<k, and call this equation (2).As before, we multiply equations (1) and (2) together, giving
(a^{2}+b^{2}+c^{2}+d^{2})(A^{2}+B^{2}+C^{2}+D^{2}) = mk^{2}p, and from the Quaternion Identity, we get
(aA+bB+cC+dD)^{2}+(aBbA+cDdC)^{2}+(aCcA+dBbD)^{2}+(aDdA+bCcB)^{2} = mk^{2}p.
As before, each of these terms is a multiple of k
Now our four squares are ((aA+bB+cC+dD)/k)^{2}, ((aBbA+cDdC)/k)^{2}, ((aCcA+dBbD)/k)^{2}, and ((aDdA+bCcB)/k)^{2}.If k is even...
The algorithm, above, fails if A=B=C=D=k/2, because then equation (2) has a righthandside of k^{2}, so m is not less than k, and it will just go on forever. Otherwise, the algorithm, above, could be used for even k as well. In any case, the following method is a little better, and it works whenever k is even.
Notice that exactly zero, two, or four of a, b, c, d are odd. Permute the numbers so that a and b are either both even or both odd, and c and d are likewise either both even or both odd. Then
((a+b)/2)^{2} + ((ab)/2)^{2} + ((c+d)/2)^{2} + ((cd)/2)^{2} = (k/2)pExample: find two numbers whose squares add up to 31
Start looking for an a such that (1a^{2})^{16} = 1a^{2} (mod 31)
When a=1, 1a^{2}=2, and (1a^{2})^{16}=2. No good.
When a=2, 1a^{2}=5, and (1a^{2})^{16}=5. No good.
When a=3, 1a^{2}=10, and (1a^{2})^{16}=10. No good.
When a=4, 1a^{2}=14, and (1a^{2})^{16}=14. Ahah!Now we let b^{2} = 1a^{2} = (1a^{2})^{16} (mod 31), so b = (1a^{2})^{8}.
So we start with the absolute least residuals a=4, b=(14^{2})^{8}=13, c=1, and d=0.
4^{2} + 13^{2} + 1^{2} + 0^{2} = 186 = kp = 6*31.
Since k is even, we use the second rule, and permute the numbers so the first two and the last two have the same parity:
4, 0, 13, and 1.
Now calculate the sums and differences, divided by 2:
(4+0)/2, (40)/2, (13+1)/2, and (131)/2 are 2, 2, 7, and 6.
2^{2} + 2^{2} + 7^{2} + 6^{2} = 93 = kp = 3*31.Since k is odd, we use the first rule, and calculate the absolutely least residuals of
2, 2, 7, 6 mod 3  they are 1, 1, 1, 0.
Now, from the quaternion identity, find
(aA+bB+cC+dD)/k = (2*1+2*1+7*1+6*0)/3 = 1
(aBbA+cDdC)/k = (2*12*1+7*06*1)/3 = 2
(aCcA+dBbD)/k = (2*17*1+6*12*0)/3 = 1
(aDdA+bCcB)/k = (2*06*1+2*17*1)/3 = 5
1^{2} + 2^{2} + 1^{2} + 5^{2} = kp = 31, so k=1, and we're done!
Existence
if p is a prime of the form 4k+1, then there is a number a such that a^{2} = 1 (mod p).
Proof:
Let h, which stands for "half", be given by h=(p1)/2. E.g. if p=13 then h=6.
Let a be any number not divisible by p, such that a^{h}=1. (From the "Polynomial Roots" discussion, below, we see that half of all possible values of a between 0 and p have this property.)
Let b=a^{h/2}, in which case b^{2}=1 (mod p).. . . . . . reference Legendre Symbol
Method
To find a number, b, such that b^{2} = 1 (mod p), follow this procedure:
Test possible values of a to find one such that a^{h}=1.
Then let b=a^{h/2}.
In practice, the smallest value of a with this property will be a prime, so you can save a little time by making a list of primes, and checking them first.
Consider the polynomial x^{p1}  1 = 0 (mod p).
By Fermat's Little Theorem, this polynomial has p1 roots, namely {1, 2, 3, ..., p1}.
Let h, which stands for "half", be given by h=(p1)/2. E.g. if p=13 then h=6.
Observe that x^{p1}1 can be factored as (x^{h}1)(x^{h}+1).
The polynomial x^{h}1=0 (mod p) can have at most h roots, and the other factor of x^{p1}1 tells us that
the polynomial x^{h}+1=0 (mod p) can have at most h roots.For any a, a^{h}=1 or a^{h}=1. In fact, for exactly half of the possible choices of a, a^{h}=1, and for the other possible choices of a, a^{h}=1.
a^{p1} = 1 (mod p), where p is prime and a≠0 (mod p)
Proof: take a to be any number 0 < a < p.
Consider the set {a, 2a, 3a, ..., (p1)a}.
No two elements of this set are congruent mod p, so their residues mod p must be {1, 2, 3, ..., p1} in some order.
By taking the products of the numbers in the two sets,
a^{p1}(p1)! = (p1)! (mod p)
Since (p1)! (mod p) is not zero, we can divide both sides by (p1)!, proving Fermat's Little Theorem.
First, let's look at the left side of this equation.
Let a and b be such that 0 <= a <= (p1)/2, and 0 <= b <= (p1)/2.
If a^{2}b^{2}=0 (mod p) then (ab)(a+b)=0 (mod p). It's not possible for a+b=0 (mod p) so a=b.
So the squares of {0, 1, ..., (p1)/2} have (p+1)/2 distinct residues.Now, for the right side of the equation.
By the similar reasoning, the numbers {10^{2}, 11^{2}, ..., 1((p1)/2)^{2}} have (p+1)/2 distinct residues.Since altogether there can be no more than p distinct residues, the two sets of residues overlap, and there must be one in common between the two sets.
So called because if u=abicjdk and U=A+Bi+Cj+Dk, uU = uU. Squaring both sides gives:
(a^{2}+b^{2}+c^{2}+d^{2})(A^{2}+B^{2}+C^{2}+D^{2}) = (aA+bB+cC+dD)^{2}+(aBbA+cDdC)^{2}+(aCcA+dBbD)^{2}+(aDdA+bCcB)^{2}
Thus, the product of any two numbers, each of which can be expressed as the sum of four squares, can also be expressed as the sum of four squares.
The product of any two numbers, each of which can be expressed as the sum of two squares, can also be expressed as the sum of two squares:
(a^{2}+b^{2})(A^{2}+B^{2}) = (aA+bB)^{2}+(aBbA)^{2}
A variant of the Complex Product Identity shows that the product of two numbers, each of which has the form a^{2}+3b^{2} can also be expressed in the form a^{2}+3b^{2}.
(a^{2}+3b^{2})(A^{2}+3B^{2})=(aA3bB)^{2}+3(aB+bA)^{2}
This identity is used in the proof of the n=3 case of Fermat's Last Theorem.
Dario Alpern's Four Squares page, which provided most of the information on this page, and in addition hinted that there is a simple proof that 1 isn't a quadratic residue (mod p) when p=4k+3.
Mathpages: The Jewel of Arithmetic: Quadratic Reciprocity provides an answer to this simple proof. Euler's Criterion for quadratic residue is that
a is a quadratic residue (mod p) iff a^{(p1)/2} = 1 (mod p).
Using the Legendre Symbol (a/p) to denote a^{(p1)/2} = 1 (mod p),
 (a*b/P) = (a/P)*(b/P).
 (1/P)=1.
 (1/P) = (1)^{(p1)/2}.
From this, there is a corollary:
 (1/P) = 1 if p mod 4 = 1.
 (1/P) = 1 if p mod 4 = 3.
Legendre Symbol  (a/p) = 1 iff a is a quadratic residue, i.e. nonzero square, (mod p), where p is an odd prime
Quadratic Residues  Modular Arithmetic and the Quadratic Equations ax^2+bx+c=0, a(p1)/2 ≡�1 (mod P), a is a quadratic residue (mod p) iff a^((p1)/2)≡1 (mod p), Legendre Symbol, Prime of form 4k+1 is a "Pythagorean Prime" because it can be expressed as the sum of two squares, Legendre (2/p) = (p^21)/8 (mod p) attributed to Gauss, quadratic reciprocity law, computing sqrt (mod p) using TonelliShanks algorithm
Irrationality Proofs  Proofs that π and e are irrational.
Prove that the area of a right triangle with integer sides is not a perfect square.
Arithmetic Sequence of Perfect Squares, page 3  If a^{2}, b^{2}, c^{2} are in arithmetic sequence, why is their constant difference a multiple of 24? Look at the second answer to this question for the relationship between Pythagorean Triples and this arithmetic sequence of squares.
Fermat's Theorems, and, in particular, the n=3 case of Fermat's Last Theorem
Octonion: contains a description of Degen's Eight Square Identity, which is the 8tuple analog to the quaternion and complex product identities described in this page.
Triangle Area Puzzle  What is the area of triangle whose sides are sqrt(61), sqrt(153), sqrt(388)? 61, 153, 388 are each sum of two squares, so the sides of this triangle can be drawn as the hypotenuses of three right triangles.
The webmaster and author of this Math Help site is Graeme McRae.