Express Integer as Sum of Four Squares This is kind of a follow-up to the question I posted here about expressing integers as the sum of two squares. Is there a similar general method for expressing integers as the sum of four squares? I believe the Lagrange's Four-Square Theorem states that all positive integers are expressible as the sum of four squares of integers, but how do you find these numbers. As an example consider the value $1638$. How can we find the four squares?
Similar to the Brahmagupta-Fibonacci two-square identity. Euler has a four square identity which involves the sum of 4 squares:

$$(a_1^2+a_2^2+a_3^2+a_4^2)(b_1^2+b_2^2+b_3^2+b_4^2) =\\ 
\quad(a_1b_1 - a_2b_2 - a_3b_3 - a_4b_4)^2 + (a_1b_2+a_2b_1+a_3b_4-a_4b_3)^2
+(a_1b_3 - a_2b_4 + a_3b_1 + a_4b_2)^2 + (a_1b_4 + a_2b_3 - a_3b_2 + a_4b_1)^2$$

Factor $1638$ as products of any small factors you know how to represent as sum of 4 squares. Repeat apply the formula will allow you to represent $1638$ itself as sum of 4 squares.

For example, let's say we have factored $1638$ as $2\cdot 3^2 \cdot 7 \cdot 13$, we have:

$$\begin{align}
  &amp; 2\cdot 3^2 \cdot 7 \cdot 13\\ 
= &amp; (1^2+1^2+0^2+0^2)(1^2+1^2+1^2+0^2)^2(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= &amp; (0^2 + 2^2 + 1^2 + 1^2)(1^2+1^2+1^2+0^2)(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= &amp; ((-3)^2 + 1^2 + 2^2 + 2^2)(2^2+1^2+1^2+1^2)(3^2+2^2+0^2+0^2)\\
= &amp; ((-11)^2+(-1)^2+2^2 + 0^2)(3^2+2^2+0^2+0^2)\\
= &amp; ((-31)^2 + (-25)^2 + 6^2 + 4^2)\\
\end{align}$$

This give you a non-trivial representation of $1638$ as $31^2 + 25^2 + 6^2 + 4^2$.

In general, there are many representations of a number as a sum of 4 squares.
There is a theorem:


  The total number of representations of a positive integer $n$ as the
  sum of four squares, representations that differ only in order and sign
  being counted as distinct, is eight times the sum of the divisors of
  $n$ that are not multiple of $4$.


The above representation is only $1$ out of $8 \sum_{d\mid 1638, 4 \nmid d} d = 34944$
ways of representing $1638$ as sum of 4 squares.
This formula seems wrong. Try letting $a_1=b_1$, $a_2=b_2$, $a_3=b_3$, $a_4=a_4$. The result should be $(a_1^2+a_2^2+a_3^2+a_4^2)$.

Euler's four-square identity says that the product of two numbers, each of which is a sum of four squares, is itself a sum of four squares. Specifically:

$$(a_1^2+a_2^2+a_3^2+a_4^2)(b_1^2+b_2^2+b_3^2+b_4^2)=
(a_1 b_1 + a_2 b_2 + a_3 b_3 + a_4 b_4)^2 +
(a_1 b_2 - a_2 b_1 + a_3 b_4 - a_4 b_3)^2 +
(a_1 b_3 - a_2 b_4 - a_3 b_1 + a_4 b_2)^2 +
(a_1 b_4 + a_2 b_3 - a_3 b_2 - a_4 b_1)^2$$