I have, for my own amusement, written the following code. However

I am in grave doubts as to the accuracy of the code in the function

"variations", having spent 50 odd years forgetting everything about

these calculations. Criticism and suggestions welcomed.

Note that the solutions presented to variations are all ordered,

but may have trailing zeroes.

/* find integral solutions to x*x + y*y + z*z ... = r*r */

/* Public domain, by C.B. Falconer, 2005-03-25 */

#include <stdio.h>

#include <limits.h>

#include <stdlib.h>

/* 181 allows operation with 16 bit ints */

#define MAXRADIUS 181

#define MAXDIMS 10

#if (INT_MAX / MAXRADIUS < MAXRADIUS)

# error "MAXRADIUS too big"

#endif

static unsigned int squares[MAXRADIUS+1];

struct soln {

int dims;

int radius;

int id;

int sv[MAXDIMS+1]; /* solution vector */

};

/* -------------- */

static void help(void)

{

printf("Usage: diophan radius [dimensions]\n"

" with radius <= %d, dimensions in 1..%d\n",

MAXRADIUS, MAXDIMS

);

exit(0);

} /* help */

/* -------------- */

/* Problem - how to compute the independant solutions

generated by negating points. 0 points can't be

negated, and two identical points are trivially

interchangeable, all of which reduce the soln count.

The negations generate reflections about some axis.

I don't remember any combinatorial maths.

*/

/* NULL input returns (and resets) accumulated sum */

static long int variations(struct soln *s)

{

long int vars;

static long int vartotal;

int d, f;

vars = vartotal;

if (!s) vartotal = 0;

else {

/* How many variations of s are available by reflections */

/* useful test group for this are (rad, dims) =

9, 3; 25, 3; 45, 3; 10, 4; 3, 9; 3, 10 */

vars = 1;

s->sv[0] = 0;

for (d = 1, f = s->dims; d <= s->dims; d++, f--) {

/* errors here on duplicated entries ?? */

if (s->sv[d]) {

/* Two possible values here, negatable */

if (s->sv[d] != s->sv[d-1]) vars *= 2;

vars *= f;

}

}

vartotal += vars;

}

return vars;

} /* variations */

/* -------------- */

static void dumpsolution(struct soln *s)

{

int d;

if (s) { /* else just accumulated total vars */

s->id++;

printf("%5d:", s->id);

for (d = 1; d <= s->dims; d++) {

printf(" %3d", s->sv[d]);

}

}

printf(" [%ld]\n", variations(s));

} /* dumpsolution */

/* -------------- */

/* brute force exhaustive search */

/* I have rediscovered the backtracking algorithm! :-) */

static void solve(int r, int dims)

{

int d; /* dimension index */

int needed[MAXDIMS+1]; /* solution criterion */

struct soln s; /* a solution develops here */

s.dims = dims; s.radius = r; s.id = 0;

for (d = 1; d <= dims; d++) s.sv[d] = 0;

needed[0] = squares[r];

s.sv[0] = 1; /* just to start the sequence */

for (d = 1; ; d++) {

for (s.sv[d] = s.sv[d-1]; s.sv[d] <= r; s.sv[d]++) {

needed[d] = needed[d-1] - squares[s.sv[d]];

if (needed[d] <= 0) { /* backtrack */

if (0 == needed[d]) dumpsolution(&s);

s.sv[d] = 0;

d--; /* use prev. dim */

if (d) continue;

return; /* all done */

}

else if (d < dims) break; /* check next dimension */

}

}

} /* solve */

/* -------------- */

int main(int argc, char **argv)

{

size_t temp;

char *errp;

unsigned int r, dims;

r = 0; dims = 2;

if ((argc < 2) || (argc > 3)) help();

if (argc > 1) {

temp = strtoul(argv[1], &errp, 10);

if ((temp > MAXRADIUS) || (errp == argv[1])

|| (temp < 1)) help();

else r = temp;

}

if (argc > 2) {

temp = strtoul(argv[2], &errp, 10);

if ((temp > MAXDIMS) || (errp == argv[2])

|| (temp < 1)) help();

else dims = temp;

}

for (temp = 0; temp <= MAXRADIUS; temp++) {

squares[temp] = temp * temp;

}

printf("solve %d %d\n", r, dims);

solve(r, dims);

dumpsolution(NULL); /* the accumulated total */

return 0;

} /* diophan main */

--

"If you want to post a followup via groups.google.com, don't use

the broken "Reply" link at the bottom of the article. Click on

"show options" at the top of the article, then click on the

"Reply" at the bottom of the article headers." - Keith Thompson