nt.number theory - Quadratic Diophantine equations solver

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