Euler's four-square product identity

Expanding the left-hand side term by term gives the sixteen products $x_i^2{y}_j^2$ with $1\le i\le 4$ and $1\le j\le 4$, each with coefficient $1$.

Squaring $z_1$ gives $z_1^2=x_1^2{y}_1^2+x_2^2{y}_2^2+x_3^2{y}_3^2+x_4^2{y}_4^2+2(x_1{x}_2{y}_1{y}_2+x_1{x}_3{y}_1{y}_3+x_1{x}_4{y}_1{y}_4+x_2{x}_3{y}_2{y}_3+x_2{x}_4{y}_2{y}_4+x_3{x}_4{y}_3{y}_4)$.

Squaring $z_2$ gives $z_2^2=x_1^2{y}_2^2+x_2^2{y}_1^2+x_3^2{y}_4^2+x_4^2{y}_3^2+2(-x_1{x}_2{y}_1{y}_2-x_1{x}_3{y}_2{y}_4+x_1{x}_4{y}_2{y}_3+x_2{x}_3{y}_1{y}_4-x_2{x}_4{y}_1{y}_3-x_3{x}_4{y}_3{y}_4)$.

Squaring $z_3$ gives $z_3^2=x_1^2{y}_3^2+x_2^2{y}_4^2+x_3^2{y}_1^2+x_4^2{y}_2^2+2(x_1{x}_2{y}_3{y}_4-x_1{x}_3{y}_1{y}_3-x_1{x}_4{y}_2{y}_3-x_2{x}_3{y}_1{y}_4-x_2{x}_4{y}_2{y}_4+x_3{x}_4{y}_1{y}_2)$.

Squaring $z_4$ gives $z_4^2=x_1^2{y}_4^2+x_2^2{y}_3^2+x_3^2{y}_2^2+x_4^2{y}_1^2+2(-x_1{x}_2{y}_3{y}_4+x_1{x}_3{y}_2{y}_4-x_1{x}_4{y}_1{y}_4-x_2{x}_3{y}_2{y}_3+x_2{x}_4{y}_1{y}_3-x_3{x}_4{y}_1{y}_2)$.

The terms $x_i^2{y}_j^2$ occurring in steps 1.2 to 1.5 are the sixteen products of step 1.1, each occurring once: $z_1^2$ contributes the pairs $(i,j)$ with $i=j$, $z_2^2$ those with $\{i,j\}=\{1,2\}$ or $\{3,4\}$, $z_3^2$ those with $\{i,j\}=\{1,3\}$ or $\{2,4\}$, and $z_4^2$ those with $\{i,j\}=\{1,4\}$ or $\{2,3\}$.

The remaining terms of steps 1.2 to 1.5 cancel in six pairs of coordinate pairs: $2x_1{x}_2{y}_1{y}_2$ and $2x_3{x}_4{y}_3{y}_4$ from $z_1^2$ against their negatives in $z_2^2$; $2x_1{x}_3{y}_1{y}_3$ and $2x_2{x}_4{y}_2{y}_4$ from $z_1^2$ against their negatives in $z_3^2$; $2x_1{x}_4{y}_1{y}_4$ and $2x_2{x}_3{y}_2{y}_3$ from $z_1^2$ against their negatives in $z_4^2$; $2x_1{x}_4{y}_2{y}_3$ and $2x_2{x}_3{y}_1{y}_4$ from $z_2^2$ against their negatives in $z_3^2$; $2x_1{x}_3{y}_2{y}_4$ and $2x_2{x}_4{y}_1{y}_3$ from $z_4^2$ against their negatives in $z_2^2$; and $2x_1{x}_2{y}_3{y}_4$ and $2x_3{x}_4{y}_1{y}_2$ from $z_3^2$ against their negatives in $z_4^2$.

Adding steps 1.2 to 1.5 and using steps 2.1 and 2.2, the sum $z_1^2+z_2^2+z_3^2+z_4^2$ equals the sixteen products of step 1.1, which is the left-hand side; since the computation used only the ring axioms, it is an identity of polynomials with integer coefficients and holds for every choice of the eight integers, negative or zero included.