I'll just give the following results :

- Every number is the sum of at most 4 non null squares (Lagrange theorem)
- Every number which is not of the form 4
^{n}(8k+7) is the sum of at most 3 non null squares - 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

- 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)² - 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. - 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... - Let N = 2
^{t}∏p^{m}∏q^{2n}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. - 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 :

A_{0}= 1, B_{0}= 0, C_{0}= -p, N = [√p]

Repeat until A_{n}+ C_{n}= 0 (stop is guaranted only if p = 4k + 1 is a prime) :

m_{i}= [(|B_{i}| + N)/|A_{i}|], A_{i+1}= A_{i}m_{i}² + 2B_{i}m_{i}+ C_{i}, B_{i+1}= A_{i}m_{i}+ B_{i}, C_{i+1}= A_{i}

Then p = |A_{n}|² + |B_{n}|².

For instance p = 269 : N = 16, A_{0}= 1, B_{0}= 0, C_{0}= -269

m_{0}= 16, A_{1}= -13, B_{1}= 16, C_{1}= 1

m_{1}= 2, A_{2}= 13, B_{2}= -10, C_{2}= -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²

- To get all decompositions as sum of two squares of
N = 2
^{t}∏p^{m}∏q^{2n}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}∏q^{n}with e = {± i, ± 1}, and for each factor 0 ≤ r ≤ mThen 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 = 5

^{4}×13^{4}×17^{4}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 NA Javascript to give all decompositions of N with this method.

- 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 = 5^{a}×13^{b}×17^{c}×29^{d}... 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 = 5^{18}

- P = 20 = 1×20 = 2×10 = 4×5 = 2×2×5 give : 5^{19}, 5^{9}×13, 5^{4}×13^{3}, 5^{4}×13×17.

The smallest is 5^{4}×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² = 5^{20}and L² = 5^{6}×13^{2}.

L² = 5^{6}× 13^{2}is the smallest and gives hypothenuse L = 5^{3}×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 R^{2}= 5^{a}×13^{b}×17^{c}... 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 5^{a}×13^{b}..., 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² = 5^{5}= 3125, 24 representations

R² = 5^{4}×13 = 8125, 40 representations

R² = 5^{2}×13^{2}= 4225, 36 representations

R² = 5^{2}×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.