/* jiw 

Re: Subject: a^2+b^2 = c^2+d^2
       Date: 17 Jul 2006 22:54:52 -0700
       From: titus_piezas@yahoo.com
 Newsgroups: sci.math.num-analysis, sci.math

Counting the number of positive-integer a,b,c,d tuples, with a,b
coprime and c,d coprime, a # c, such that a^2+b^2 = c^2+d^2 < k^2,
for various k, or in this program version, a^3+b^3 = c^3+d^3 < k^3.

Copyright 2006 James Waldby.  Offered without warranty
under GPL terms as at http://www.gnu.org/licenses/gpl.html

*/
#define zz fprintf(stderr,"%s@%3d \n",__FILE__,__LINE__);
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
//================================================
#define HITE 3
#if HITE==2
#define POWER(x) x*x
#define ROOT(x) sqrt(x)
#endif
#if HITE==3
#define POWER(x) x*x*x
#define ROOT(x) cbrt(x)
#endif
#if HITE==4
#define POWER(x) x*x*x*x
#define ROOT(x) sqrt(sqrt(x))
#endif
#define KMULTI 1
#define KDELT 20
#define KONE 20
#define ZMAX 40000001
#define ZHAT ((7*ZMAX)/5)
#define DEB 0
//================================================
// gcd(a,b) by Euclid's algorithm
int gcd(int a, int b) {
  int r = a%b;
  while (r) {
    a=b; b=r; r=a%b;
  }
  return b;
}
//================================================
int main(int argc, char *argv[]) {

  int a, b, c, d, g, gab, h;
  int i, ii, j, k, t, u, vnext, z;
  int *cl, *cv;

  if(DEB)printf ("\nbegin\n");
  printf ("\npower(2)=%d,  root(4096)=%f\n", POWER(2),ROOT(4096));
  cl = malloc(ZHAT*sizeof(int));
  cv = malloc(ZHAT*sizeof(int));
  memset (cl, 0, ZHAT*sizeof(int)); // Zero the links array
  vnext = ZMAX;
  for (i=1; ; ++i) {
    ii = POWER(i);
    if (ii > ZMAX) break;
    for (j=i; ; ++j) {
      t = ii + POWER(j);
      if (t > ZMAX) break;
      if (cl[t]) {		// Add [i j] to list t
	cv[vnext] = i;
	cl[vnext] = cl[t];
	cl[t] = vnext;
	++vnext;
	if (vnext >ZHAT) {
	  printf ("\nOverflowing arrays, vnext = %d\n", vnext);
	  exit (0);
	}
      } else {			// Begin list t with [i j]
	cl[t] = t;
	cv[t] = i;
      };
    }
  }
  // For various z, count cases where a^n+b^n=c^n+d^n with (a,b,c,d)=1
  h = 1;
  j = 0;
  k=KONE;
  z=POWER(k);
  do {
    for (; h<z; ++h) {
      if (cl[h]==0 || cl[h]==h) // Check if count is 0 or 1
	continue;
      t = h;
      if(DEB)printf ("\nAt h=%4d  cl_h %d   cv_h %d   j=%d   z=%d\n",
	      h, cl[h], cv[h], j, z);
      do {			// Process list h
	if(DEB)printf ("At h=%4d  t=%4d  j=%d\n", h, t, j);
	a = cv[t];
	b = 0.5+ROOT(h-POWER(a));
	gab = gcd(a,b);

	// Inner loop: Test element t vs. following elements
	for (u=cl[t]; u != h; u=cl[u]) {
	  if(DEB)printf ("At h=%4d  t=%4d  u=%4d   j=%d\n", h, t, u, j);
	  c = cv[u];
	  d = 0.5+ROOT(h-POWER(c));
	  g = gcd(gcd(c,d), gab);
	  if(DEB)printf ("At h=%4d  abcd=%4d %4d %4d %4d  g=%d\n",
		  h, a,b,c,d, g);
	  if (g==1) {
	    ++j;
	    if(DEB)printf ("At h=%4d  abcd=%4d %4d %4d %4d  j=%d\n",
		    h, a,b,c,d, j);
	  }
	} // end of inner loop

	t = cl[t];		// Go on to next element of list
      } while (t != h);		// Until end of list
    }
    // Report total # pairs j, for "radius" k
    printf ("%4d %6d  %6.3f\n", k, j, (double)j/k);
    k = KMULTI*k + KDELT;	// Compute new "radius"
    z = POWER(k);
  } while (z < ZMAX);
  if(DEB)printf ("\ndone\n");
}