Is there software that helps list small solutions of the Diophantine equation $$x_0^2=1+x_1^2+x_2^2+\cdots+ x_n^2$$ where "small" is negotiable, but e.g. we could fix $x_0$ and and ask for the list of all solutions $x_1, \dots, x_n\in\mathbb Z$?

-
This is $$(x_0-x_1)(x_0+x_1)=x_0^2-x_1^2=1+x_2^2+\cdots +x_n^2,$$ so sort of reduces to factorization. – Robin Chapman Jun 10 '10 at 20:17
@Will, I wish to distinguish $x_j$ and $-x_j$. Please do not bother to write anything in C++; I just want to know whether there is anything already available as part of Mathematica ot other similar package. @Robin, thanks, reducing factorization is helpful idea. – Igor Belegradek Jun 10 '10 at 21:01
For $n=2$ you can try alpertron.com.ar/QUAD.HTM for various values of $x_0$. – Eric Rowell Jun 10 '10 at 21:12
SquaresR[d, n] gives the number of ways r_d (n) to represent the integer n as a sum of d squares. PowersRepresentations[n, k, p] gives the distinct representations of the integer n as a sum of k non-negative p\^th integer powers. EllipticTheta[a,u,q] gives the theta function Subscript[[CurlyTheta], a](u,q) (a=1,[Ellipsis],4). – Will Jagy Jun 10 '10 at 21:19
note the above are Mathematica commands – Will Jagy Jun 10 '10 at 21:52

Have you had a look on this tutorial: http://reference.wolfram.com/mathematica/tutorial/DiophantineReduce.html ? To take different values n < n_max into account, a simple loop could work.

-

This function in Mathematica find them all:

f[n_] := Reduce[ Total[Table[x[i]^2, {i, 1, n}]] + 1 == x[0]^2 &&
(Table[x[i], {i, 1, 3}] /. List -> LessEqual),
Table[x[i], {i, 0, n}], Integers]


Invoke with:

f[1], f[3] ... etc

-