#include <math.h>

int main(int argc, char **argv);
void solve(double a,double b,double c,double d,double e,double f);

double bb,dd;

int main(int argc, char **argv) {
  int p, m, A, C;
  int alpha, beta, E;
  int B, D;
  int RHS1, RHS2;
  double b, d;
  
  for (p=-100; p<101; p++) {
    for (m=-100; m<101; m++)
    for (A=0; A<101; A++)
    for (C=0; C<101; C++) {
      /* CASE 1: choose the top sign */
      alpha = p+m;
      beta  = p-2*m;
      E     = 3*p-(A+C);
      
      RHS1  = 3*(m*m-p*p) + A*C + A*E + C*E;
      RHS2  = A*C*E - alpha*alpha*beta;
      
      /* we don't want one or both of B and D to be zero */
      /*   NB: RHS1 is the same in cases 1 and 2         */
      if (RHS1 <= 1)
	continue;
      
      if (E!=A) {
	solve(1.0,1.0,
	      ((double)E),((double)A),
	      ((double)RHS1),((double)RHS2));
	if ((bb>0.0) && (dd>0.0)) {
	  b = round(sqrt(bb));
	  d = round(sqrt(dd));
	  if ( (abs(b*b-bb) < 0.1) && (abs(d*d-dd) < 0.1) ) {
	    printf("Found one!\n");
	    printf("p = %d\nm = %d\nA = %d\nC = %d\n",p,m,A,C);
	    printf("E = %d\nb = %g\nd = %g\n\n",E,b,d);
	  }
	}
      } else {
	if ((A!=0) && (A*RHS1==RHS2)) {
	  printf("Further checking needed...\n");
	  printf("RHS1 = %d\np = %d\nm = %d\nA = %d\nC = %d\n\n",RHS1,p,m,A,C);
	}
	if ((A==0) && (RHS2==0)) {
	  printf("Further checking needed...\n");
	  printf("RHS1 = %d\np = %d\nm = %d\nA = %d\nC = %d\n\n",RHS1,p,m,A,C);
	}
      }
      
      /* CASE 2: choose the bottom sign */
      alpha = p-m;
      beta  = p+2*m;
      E     = 3*p-(A+C);
      
      RHS1  = 3*(m*m-p*p) + A*C + A*E + C*E;
      RHS2  = A*C*E - alpha*alpha*beta;
      
      /* we don't want one or both of B and D to be zero */
      /*    NB: not really necessary, but we'll leave it */
      if (RHS1 <= 1)
	continue;
      
      if (E!=A) {
	solve(1.0,1.0,
	      ((double)E),((double)A),
	      ((double)RHS1),((double)RHS2));
	if ((bb>0.0) && (dd>0.0)) {
	  b = round(sqrt(bb));
	  d = round(sqrt(dd));
	  if ( (abs(b*b-bb) < 0.1) && (abs(d*d-dd) < 0.1) ) {
	    printf("Found one!\n");
	    printf("p = %d\nm = %d\nA = %d\nC = %d\n",p,m,A,C);
	    printf("E = %d\nb = %g\nd = %g\n\n",E,b,d);
	  }
	}
      } else {
	if ((A!=0) && (A*RHS1==RHS2)) {
	  printf("Further checking needed...\n");
	  printf("RHS1 = %d\np = %d\nm = %d\nA = %d\nC = %d\n\n",RHS1,p,m,A,C);
	}
	if ((A==0) && (RHS2==0)) {
	  printf("Further checking needed...\n");
	  printf("RHS1 = %d\np = %d\nm = %d\nA = %d\nC = %d\n\n",RHS1,p,m,A,C);
	}
      }
    }
    printf("p=%d done\n\n",p);
  }
  
  return 0;
}


void solve(double a,double b,double c,double d,double e,double f) {
  /* just be lazy and use cramer's rule */
  /* this assumes our matrix is non-singular */
  
  double detm = a*d - b*c;
  double detx1 = e*d - b*f;
  double detx2 = a*f - e*c;
  
  bb = detx1/detm;
  dd = detx2/detm;
}

