Pythagorean 5-tuples
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 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)\, .$
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).