# Sum of squares

Look at the work of Fermat, Euler, Gauss etc for details and proofs.
I'll just give the following results :
1. Every number is the sum of at most 4 non null squares (Lagrange theorem)

2. Every number which is not of the form 4n(8k+7) is the sum of at most 3 non null squares

3. The only numbers which are sum of two squares are (Fermat-Girard theorem) :
• number 2 = 1² + 1²
• The prime numbers of the form 4k + 1
• The square of prime numbers of the form 4k + 3
•  and the products of theses numbers

4. However if two numbers are each the sum of two squares, their product is also sum of two squares, in two ways :
(a² + b²)(c² + d²) = (ac + bd)² + (ad - bc)² = (ac - bd)² + (ad + bc)²

5. But a product of two sums of 3 squares is not always a sum of 3 squares :
26 = 4² + 3² + 1² et 38 = 5² + 3² + 2² but 26×36 = 988 = 4×(8×30 + 7) is of the form (2) and can't be sum of 3 squares.

6. Mention also the Euler relation for sums of 4 squares :
(a² + b² + c² + d²)(x² + y² + z² + t²) = (ax + by + cz + dt)² + (ay - bx + ct - dz)² + (az - cx + dy - bt)² + (at - dx + bz - cy)²
and all analogue formulaes obtained by changing b to -b etc...

7. Let N = 2t∏pm∏q2n the prime decomposition of N, where we separate factor 2, prime factors p of the form 4k + 1 and prime factors q of the form 4k + 3.
N is sum of two squares if and only if the exponents of the 4k+3 primes are even
Then the number of ways to write N as sum of two squares is r(N) = 4×∏(m + 1) and depends only on the exponents of the 4k + 1 primes.
The representations (±x, ±y) and (±y, ±x) count as 8 different ways, except when x = 0, y = 0 or x = y.
The number of really different decompositions of N = x² + y², x ≥ y ≥ 0 is :

 d(N) = int[(1 + ∏(m + 1))/2]

Example :
Number of representations of 1521000 = 8×3²×5³×13² is 4×(3 + 1)(2 + 1) = 48 (only 5 and 13 are of the form 4k + 1, with exponents 3 and 2). The number of different decompositions is [(1 + (3 + 1)(2 + 1))/2] = 6.

8. To get the decomposition as sum of 2 squares of a prime 4k + 1, we can do as follows :

• By trying x = 0, 1, 2... until we find a square p - x². obvious, simple, not very efficient.

• Or by repeating :
y = [_ √(p - x²) _] and x = [¯ √(p - y²) ¯] starting from x = 0 until a null remainder or x > y.
With [_ z _] largest integer n ≤ z and [¯ z ¯] smallest integer n ≥ z.
For instance p = 269 x = 0, y = 16, x = 4, y = 15, x = 7, y = 14, x = 9, y = 13, x = 10 finished : 269 = 10² + 13²

• Or also with the Legendre method :
A0 = 1, B0 = 0, C0 = -p, N = [√p]
Repeat until An + Cn = 0 (stop is guaranted only if p = 4k + 1 is a prime) :
mi = [(|Bi| + N)/|Ai|], Ai+1 = Aimi² + 2Bimi + Ci, Bi+1 = Aimi + Bi, Ci+1 = Ai
Then p = |An|² + |Bn|².
For instance p = 269 : N = 16, A0 = 1, B0 = 0, C0 = -269
m0 = 16, A1 = -13, B1 = 16, C1 = 1
m1 = 2, A2 = 13, B2 = -10, C2 = -13 finished : 269 = 13² + 10²

• And finally with the Heath-Brown method :
With p = 4k + 1 let x = k, y = 1, z = 1 then repeat
if x > y + z : x = x - y - z, y = y, z = 2y + z
if x < y + z : x = y + z - x, y = x, z = 2x - z
(invariant 4xy + z² = p, x - y + z > 0)
until x = y. Then p = 4x² + z²
For instance with p = 29 = 4*7 + 1 x = 7, y = 1, z = 1, y + z = 2
x = 5, y = 1, z = 3, y + z = 4
x = 1, y = 1, z = 5 finished 29 = 4 + 25 = 2² + 5²

9. To get all decompositions as sum of two squares of N = 2t∏pm∏q2n with p = u² + v² being the primes of the form 4k + 1 and q the primes of the form 4k + 3, let's define the complex number :

A + iB = e(1 + i)t ∏(u + iv)r(u - iv)m-r ∏qn with e = {± i, ± 1}, and for each factor 0 ≤ r ≤ m

Then N = A² + B² is a way of writing N as sum of two 2 squares. We get all of them by varying e and each r. We get only the decompositions (non equivalent) by fixing e and some constraints on the r.

Example : N = 1521000 = 8×3²×5³×13², A + iB = e(1 + i)³(2 + i)r(2-i)3-r(3 + 2i)s(3 - 2i)2-s×3
We can fix e = -i so as e(1 + i)³ = 2(1 + i), and constraint r = {0,1,2,3} and s = {0,1,2}.
But if we choose independantly r and s, we see that r = 0, s = 0 and r = 3, s = 2 result into equivallent representations -258-1206i and -1206-258i because numbers (2 - i)³(3 - 2i)² and (2 + i)³(3 + 2i)² are conjugate.

To get only the 6 decompositions we have to restrict r = {0, 1} s = {0, 1, 2}.
r = 0, s = 0 : 6(1 + i)(2 - i)³(3 - 2i)² = -258-1206i that is N = 1206² + 258²
r = 0, s = 1 : 6(1 + i)(2 - i)³(3 + 2i)(3-2i) = 1014-702i that is N = 1014² + 702²
r = 0, s = 2 : 6(1 + i)(2 - i)³(3 + 2i)² = 1038+666i that is N = 1038² + 666²
r = 1, s = 0 : 6(1 + i)(2 + i)(2 - i)²(3 - 2i)² = 810-930i that is N = 930² + 810²
r = 1, s = 1 : 6(1 + i)(2 + i)(2 - i)²(3 + 2i)(3-2i) = 1170+390i that is N = 1170² + 390²
r = 1, s = 2 : 6(1 + i)(2 + i)(2 - i)²(3 + 2i)² = 90+1230i that is N = 1230² + 90²

Generally, we restrict one of the r to half his values. Choose a prime with an odd exponent. If every p has even exponent (N is a square or the double of a square), r shall be restricted to half its values, plus one, recursively with all exponents.

For instance N = 54×134×174 A + iB = (2 + i)r(2-i)4 -r(3+2i)s(3-2i)4 -s(4+ i)t(4-i)4 -t
r = {0, 1}, s = {0, 1, 2, 3, 4}, t = {0, 1, 2, 3, 4} : restrict r, 50 decompositions
r = 2, s = {0, 1}, t = {0, 1, 2, 3, 4} : r fixed, restrict s, 10 decompositions
r = 2, s = 2, t = {0, 1} : r and s fixed, restrict t, 2 decompositions
r = 2 s = 2 t = 2 : r, s, t fixed, 1 decomposition
That is the expected 63 decompositions of N

A Javascript to give all decompositions of N with this method.

## Applications :

• Find the smallest number which is sum of two squares in 10 different ways?
The smallest number will have only prime factors p = 4k + 1 as small as possible N = 5a×13b×17c×29d... with a ≥ b ≥ c ≥ d... Exponents will be such as P = (a + 1)(b + 1)(c + 1)...= 19 or 20.
- P = 19 = 1×19 gives N = 518
- P = 20 = 1×20 = 2×10 = 4×5 = 2×2×5 give : 519, 59×13, 54×133, 54×13×17.
The smallest is 54×13×17 = 138125.

• Smallest hypothenuse common to 10 right triangles with integer sides.
The square of hypothenuse L² is a sum of squares in 10 different ways. The classical formulae above include L² + 0². If we want real triangles, we must discard this "décomposition", hence there are 11 and P = 21 or 22.
L² being a square (shure !), P is odd.
P = 21 = 1×21 = 3×7 L² = 520 and L² = 56×132.
L² = 56× 132 is the smallest and gives hypothenuse L = 53×13 = 1625.

Note : the smallest hypothenuse common to at least 10 right triangles might not be this one ! To find this hypothenuse, look at the products of odd numbers greater or equal to 21. 3×7 = 21 is the one found above but 3×3×3 = 27 resulting into L² = 5²×13²×17² gives a smaller hypothenuse = 5×13×17 = 1105, common to 13 right triangles, 14 with L² + 0².

• Circle going through the maximum number of integer points.
Now R2 = 5a×13b×17c... minimum with a ≥ b ≥ c... and d(R²) maximum, but R may not be an integer.
On squared A4 = 210×297 paper : 2R ≤ 210 = 42 squares, that is R² ≤ 441.
R² = 5²×13 = 325 : 24 representations (1,18) (6,17) (10,15) and symetrics.
R² = 5³ = 125 : 16 representations only (2,11) (5,10) and symetrics.
If R is an integer, then R itself is 5a×13b..., and not only R².
R ≤ 21 gives R = 5, R² = 25 and 12 representations only (0,5) (3,4) and symetrics.

With millimeter squared paper, 2R ≤ 210mm and R² ≤ 11025. The only posibilities are :
R² = 55 = 3125, 24 representations
R² = 54×13 = 8125, 40 representations
R² = 52×132 = 4225, 36 representations
R² = 52×13×17 = 5525, 48 representations (7,74) (14,73) (22,71) (25,70) (41,62) (50,55) and symetrics.
The maximum is that last value.