Pythagorean 5-tuples What is the solution of the equation $x^2+y^2+z^2+t^2=w^2$ in polynomials over C ("Pythagorean 5-tuples")? 
There are simple formulas describing Pythagorean n-tuples for n=3,4,6:


*

*n=3. The formula for solutions of $x^2+y^2=z^2$ [4]:


$x=d(p^2−q^2)$,
$y=2dpq$, 
$z=d(p^2+q^2),$
where p,q,d are arbitrary polynomials. 


*

*n=4. Similarly, all Pythagorean 4-tuples of polynomials are given by the identity [4, Theorem 2.2], [3, Theorem 7.1]:


$(p^2+q^2-r^2-s^2)^2+(2pr+2qs)^2+(2ps-2qr)^2=(p^2+q^2+r^2+s^2)^2$


*

*n=6. The following identity produces Pythagorean 6-tuples [3, Theorem 7.2]:


$(m^2-n^2)^2+(2mn)^2+(2(n_0m_1-m_1n_0+m_3n_2-m_2n_3))^2+$
$(2(n_0m_2-m_2n_0+m_1n_3-m_3n_1))^2+$
$(2(n_0m_3-m_3n_0+m_2n_1-m_1n_2))^2=(m^2+n^2)^2$
where $m=(m_1,m_2,m_3,m_4)$, $n=(n_1,n_2,n_3,n_4)$, and $mn$ is the usual scalar product.
These identities are somehow related to sl(2,R), sl(2,C), sl(2,H), but the case n=5 is missing in this description [3, Theorems 7.1 and 7.2].  
The above formuli describe also Pythagorean n-tuples of integers. There are another descriptions of those; see [2] and [5, Chapter 5] for n<10, [6] for n<15, [1, Theorem 1 in Chapter 3] for arbitrary n, and also the answers below. 
There are reasons to believe that Pythagorean 5-tuples cannot be described by a single polynomial identity. Thus identities producing a "large" set of solutions are also of interest, like 
$(-p^2+q^2+r^2+s^2)^2+(2pq)^2+(2pr)^2+(2ps)^2=(p^2+q^2+r^2+s^2)^2$
This identity does not give all solutions because it produces only reducible polynomials y,z,t (once p,q,r,s are nonconstant). Examples of solutions which cannot be obtained by the approaches in the answers below are also interesting.
Given a solution (x,y,z,t,w), methods to construct a new solution (x',y',z',t',w') are also of interest. For instance, 
$x'=w+y+z$,
$y'=w+z+x$,
$z'=w+x+y$,
$t'=t$,
$w'=2w+x+y+z$
[see the answer of Ken Fan below for generalizations to other n] or
$x'+iy'+jz'+kt'=q(x+iy+jz+kt)q$,
where $q$ is arbitrary polynomial with quaternionic coefficients [see the answer of Geoff Robinson below].
--
[1] L.J. Mordell, Diophantine Equations, Academic Press, London, 1969
[2] D. Cass, P.J. Arpaia, MATRIX GENERATION OF PYTHAGOREAN n-TUPLES, Proc. AMS 109:1 (1990), 1-7
[3] J. Kocik, Clifford Algebras and Euclid’s Parametrization of Pythagorean Triples,  Adv. Appl. Clifford Alg. 17:1 (2007), 71-93 
[4] R. Dietz, J. Hoschek and B. Juttler, An algebraic approach to curves and surfaces on the sphere
and on other quadrics, Computer Aided Geom. Design 10 (1993) 211-229
[5]  V. Kac, Infinite-Dimensional Lie Algebras (3rd edn. ed.), CUP, 1990
[6] E. Vinberg, The groups of units of certain quadratic forms  (Russian), Mat. Sbornik (N.S.) 87(129) (1972), 18–36 
 A: In the early 90s, I took a class from Victor Kac and in it, he explained a very beautiful way of generating all the primitive solutions for the Pythagorean equation for a sum of $n-1$ perfect squares is equal to a perfect square where $n$ can be 3, 4, 5, ..., 10.  Unfortunately, I do not know where in the literature this is described in detail.  It might be in his Infinite Dimensional Lie Algebras book, but I don't know.
The idea is to realize the solutions as the isotropic roots for a certain root system.
Consider the lattice ${\Bbb Z}^n$ with bilinear form $-x_0y_0 + x_1y_1 + \cdots + x_{n-1}y_{n-1}$ and standard basis $v_0,\ v_1,\ v_2,\ \dots,\ v_{n-1}$.
Change basis to:
$\alpha_1 = v_1 - v_2$, $\alpha_2 = v_2 - v_3$, $\dots$, $\alpha_{n-2} = v_{n-2} - v_{n-1}$, and $\alpha_{n-1} = v_{n-1}$.
If $n \geq 4$, let $\alpha_n = -v_0 - v_1 - v_2 - v_3$.
If $n=3$, let $\alpha_n = -v_0 - v_1 - v_2$.
The corresponding Cartan matrix $a_{ij} = \frac{2(\alpha_i, \alpha_j)}{(\alpha_i, \alpha_i)}$ is represented by the diagram:
    
 (source)
Then the set of primitive solutions to the equation $x_0^2 = x_1^2 + x_2^2 + \cdots + x_{n-1}^2$ is the orbit under the corresponding Coxeter group of $(1, 1, 0, \dots, 0)$ if $n < 10$.  If $n=10$, then you have to add the orbit $(3, 1, 1, 1, \dots, 1)$ to get them all.
One doesn't need the theoretical machinery to prove the result.  One can just construct the matrices and use a descent argument to show that it works.
For the Pythagorean triple case, for instance, you take the orbit of the vector $(1, 1, 0)$ under the action of the group generated by the matrices:
$$
\begin{pmatrix}
1 & 0 & 0 \\
0 & 0 & 1\\
0 & 1 & 0
\end{pmatrix},
$$
$$
\begin{pmatrix}
\pm 1 & 0 & 0 \\
0 & \pm 1 & 0 \\
0 & 0 & \pm 1
\end{pmatrix},
$$ and
$$
\begin{pmatrix}
3 & 2 & 2 \\
2 & 1 & 2 \\
2 & 2 & 1
\end{pmatrix}.
$$
For $n=4$, you can use the matrices that permute the appropriate variables, change the sign of any variable, and the following:
$$
\begin{pmatrix}
2 & 1 & 1 & 1 \\
1 & 0 & 1 & 1 \\
1 & 1 & 0 & 1 \\
1 & 1 & 1 & 0
\end{pmatrix}.
$$
I'm sorry I cannot give a reference to the literature...I hope someone else is able to.
A: The construction below shows that such tuples really hardly could be written in unique way. Maybe some exclusion from general case is due to the Hurwitz theorem about sum of squares and noncommutativity of quaternions. This construction is modification of Pietro Majer and Geoff Robinson suggestions. 
Let us consider quaternion $q = t + x i + y j + z k$ and a constant quaternion $c$. I am not sure after reading of the question, if norm should be a square of some integer or not. It may be Hurwitz quaternion and simple nontrivial example is $c=(1+i+j+k)/2$. Now let us consider product $s = q c q$ as polynomial of $t, x, y, z$. 
We have $|s| = |c| |q|^2 = |c| (t^2+x^2+y^2+z^2)^2 $, but $|s|=s_0^2+s_1^2+s_2^2+s_3^2$ for $s = s_0 + s_1 i + s_2 j +s_3 k$, where $s_0, s_1, s_2, s_3$ by definition are second order polynomial of $t, x, y, z$.
So different $c$ produces different tuples and $c=1$ produces solution mentioned by  Pietro Majer and Geoff Robinson and in Mordell book.
[edit] To produce more general solution for second order polynomials it is possible to use some modification. Let's consider three constant quaternions $a, b, c$ with modules are squares of some integers. It may be done using Kac's method or something else. Then new solution is $a q c q b$.
A: Another way to generate parametrisations is to repeatedly apply a parametrisation of triples:
Say we start with 
$$(x(s), y(s)) = (\frac{2s}{s^2 + 1}, \frac{s^2 - 1}{s^2 + 1}),$$ satisfying $$x(s)^2 + y(s)^2 = 1.$$ 
To parametrise solutions to $x^2 + y^2 + z^2 = 1$, write $x^2 + y^2 = 1 - z^2$ which we set equal to $w^2$, so that a parametrised solution $(w(t), z(t))$ to $w^2 + z^2 = 1$ gives rise to a parametrised solution $$(w(t)x(s), w(t)y(s), z(t)) = (\frac{2t}{t^2 + 1}\cdot\frac{2s}{s^2 + 1}, \frac{2t}{t^2 + 1}\cdot\frac{s^2 - 1}{s^2 + 1}, \frac{t^2 - 1}{t^2 + 1})$$ to $x^2 + y^2 + z^2 = 1$.
Repeating this, we get your case:
$$(\frac{2u}{u^2 + 1}\cdot\frac{2t}{t^2 + 1}\cdot\frac{2s}{s^2 + 1}, \frac{2u}{u^2 + 1}\cdot\frac{2t}{t^2 + 1}\cdot\frac{s^2 - 1}{s^2 + 1}, \frac{2u}{u^2 + 1}\cdot\frac{t^2 - 1}{t^2 + 1}, \frac{u^2 - 1}{u^2 + 1})$$ 
parametrises
$$x^2 + y^2 + z^2 + w^2 = 1$$
(to go to your original problem you can throw in another parameter to homogenise, clear denominators and multiply everything by another parameter).
If you apply the same to the parametrisation $(\cos(t), \sin(t))$ of the circle, you get a parametrisation of the 2-sphere by sperical coordinates.
A: A parametrization of solutions is
$x=2ad$
$y=2bd$
$z=2cd$
$t=a^2+b^2+c^2-d^2$
$w=a^2+b^2+c^2+d^2\, .$
It is easy to see that this generates all rational solutions if $a, b, c, d$ are rational numbers, and (consequently) all integers solutions up to a similarity factor, if $a, b, c, d$ are integers.
[edit] for what it's worth, here's something more symmetric (with the same rmk as before)
$x=-a^2+b^2+c^2+d^2-2a(b+c+d)$
$y=\phantom{-}a^2-b^2+c^2+d^2-2b(a+c+d)$
$z=\phantom{-}a^2+b^2-c^2+d^2-2c(a+b+d)$
$t=\phantom{-}a^2+b^2+c^2-d^2-2d(a+b+c)$
$w=\phantom{-}2(a^2+b^2+c^2+d^2)\, .$
A: One way to generate integer solutions is as follows: Let $ v = (p + qi + rj + sk)$, where $p,q,r,s$
are rational integers and $\{1,i,j,k\}$ is the usual $\mathbb{R}$-basis for the algebra of quaternions. Let $v^{\prime}= (p -qi -rj - sk)$ and 
$w = vv^{\prime} = (p^{2}+q^{2}+r^{2}+s^{2}).$ 
If we expand $v^{2}$ in the form $x + yi + zj +tk$ for integers $x,y,z,t$, then we do
have $x^{2} + y^{2} + z^{2} + t^{2} = w^{2}$. This is somewhat analogous to generating 
the Pythagorean triple $(x^{2}-y^{2})^{2}$ + $(2xy)^{2}$ = $(x^{2}+y^{2})^{2}$ by taking 
$({\rm Re}(x+iy)^{2})^{2} + ({\rm Im}(x+iy)^{2})^{2} = (x^{2}+y^{2})^{2}$. It may be
better to work with Hurwitz quaternions for this question. 
Answer extended following the answer of Alex qubeat below, and slight rephrasing of the
original question to place more emphasis on polynomials: At least in the context
of polynomials in $\mathbb{R}[u]$ (it is easy to run out of letters in this game, so $u$
denotes an indeterminate here), more solutions may be manufactured using the fact that the
solutions to $x(u)^{2} +y(u)^{2} + z(u)^{2} + t(u)^{2} = w(u)^{2}$ with 
$x(u),y(u),z(u),t(u),w(u) \in \mathbb{R}[u]$ have a semigroup structure.
Alex's answer combines a fixed solution with an essentially ``constant" solution, to produce 
other solutions, but solutions can be combined in other ways. Let $\mathbb{H}$ denote the algebra of real quaternions. Then the map $N: \mathbb{H}[u] \to \mathbb{R}[u]$ with 
$N(x(u) +iy(u)+jz(u) + kt(u)) =x(u)^{2} +y(u)^{2} + z(u)^{2} +t(u)^{2}$ (for real polynomials
$x(u),y(u),z(u),t(u))$ is multiplicative. The polynomials $p(u) \in \mathbb{H}[u]$ such 
that $N(p(u))$ is a non-zero square in $\mathbb{R}[u]$ are closed under multiplication.
As in the integral case, one way to ensure that $N(x(u) + iy(u) + jz(u) + kt(u))$ is
a square is to take $x(u) + iy(u) +jz(u)+kt(u)$ to be of the form $(a(u)+ib(u)+jc(u)+kd(u))^{2}$
for real polynomials $a(u),b(u),c(u),d(u)$, but it is also clear from this discussion
that $N(p_{1}(u)p_{2}(u) \ldots p_{n}(u))$ is a square in $\mathbb{R}[u]$ as long 
as any given $p_{i}(u)\in \mathbb{H}[u]$ occurs an even number of times in the (non-commuting) product. In fact, it is permissible to count the total number of occurences of $\mathbb{H}$-conjugates of any given $p_{i}(u)$, that is elements of the form $v^{-1}p_{i}(u)v$, where $v$ is a non-zero element of $\mathbb{H}$,  since $N$ takes a constant value on the $\mathbb{H}^{\times}$-conjugacy class of $p_{i}(u)$.

Further remarks. May 6: The polynomial ring $\mathbb{H}[u]$ has a natural division ring of fractions,
$\mathbb{H}(u)$, which is isomorphic to a certain ring of $2 \times 2$ matrices over $\mathbb{C}(u)$
(I can't get the latex right for matrices, but the ring should be clear- in particular, 
the determinants of elements in the ring are elements of $\mathbb{R}(u)$). There is a natural
ring homomorphism $\sigma$ from $\mathbb{H}(u)$ to ${\rm SO}(3,\mathbb{R}(u))$, obtained by letting
$\mathbb{H}(u) \backslash \{0\}$ act by conjugation of the $\mathbb{R}(u)$-span of $\{i,j,k\}.$
The question about polynomials amounts to determining $S\sigma$, where $S$ is the set of elements
of $\mathbb{H}(u)$ whose determinants are squares in $\mathbb{R}(u)$. This is the same as
$T\sigma$, where $T$ is the set of elements of determinant $1$, since the non-zero elements
of $\mathbb{R}(u)$ acts trivially by conjugation on $\mathbb{H}(u)$. Now it is clear
that $T\sigma$ is a normal subgroup of ${\rm SO}(3,\mathbb{R})$, and that the factor
group is an elementary Abelian $2$-group (that is, an Abelian
group of exponent $2$). I do not see at present how to calculate the size of this group,
but record this line of thinking in case anyone else can exploit it further.

More remarks, May 13: I think that $\mathbb{H}(u)$ has an involutory automorphism $\sigma$
which fixes $\mathbb{R}(u)$ elementwise, and which has $i \sigma = -i$,$j\sigma = k$ and 
$k\sigma =j.$ Then $N(p(u)\sigma) = N(p(u))$ for all $p(u) \in \mathbb{H}(u)$.
Hence, for any $p(u) \in \mathbb{H}(u)$, the element $q(u) = p(u).p(u)\sigma$
will have $N(q(u)) =w(u)^{2}$ for some $w(u) \in \mathbb{R}(u)$, so this will yield
a new type of solution (at least in the context of this answer).
A: In fact, the equation:
$X^2+Y^2+Z^2+R^2=W^2$
Solutions look like this:
$X=2psabk^2+a^2k^2s^2-ckabs^2-abk^2s^2$
$Y=2psabk^2+a^2k^2s^2-ckabs^2+2abk^2s^2$
$Z=2psabk^2+a^2k^2s^2+2ckabs^2-abk^2s^2$
$R=2p^2b^2k^2+c^2b^2s^2+b^2k^2s^2-ckb^2s^2-a^2k^2s^2+2psabk^2$
$W=2p^2b^2k^2+2psabk^2+c^2b^2s^2+b^2k^2s^2-ckb^2s^2+2a^2k^2s^2$
And the formula:
$X=2psabk^2+abk^2s^2+ckabs^2-a^2k^2s^2$
$Y=a^2k^2s^2-ckabs^2+2abk^2s^2-2psabk^2$
$Z=a^2k^2s^2+2ckabs^2-abk^2s^2-2psabk^2$
$R=2p^2b^2k^2+c^2b^2s^2+b^2k^2s^2-ckb^2s^2-a^2k^2s^2-2psabk^2$
$W=2p^2b^2k^2+c^2b^2s^2+b^2k^2s^2-ckb^2s^2+2a^2k^2s^2-2psabk^2$
Where the numbers: $p,c,s,a,k,b$ integers and set us.
Quite often, after a number of substitutions should be divided by the greatest common divisor.
