/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/27/06 (dkc)							     */
/*									     */
/*  This C program finds a and b such that (a**p + b**p)/(a + b) is a cube   */
/*  or p times a cube.	The prime factors of [(a**p+b**p)/(a+b)] must be of  */
/*  the form p**2*k+1.	p is set to 3.	Only a and b where (a**p+b**p)/(a+b) */
/*  has one distinct prime factor are output.				     */
/*									     */
/*  The output is "a, b, split".  If 2p does not divide a, b, a-b, or a+b,   */
/*  then "split" is set to 1.  If a is even, p divides a, and p**2 does not  */
/*  divide a, then an error is indicated ("error[1]" is set to 7).  b and    */
/*  a-b are treated similarly.	If a+b is even, p divides a+b, and p**3      */
/*  does not divide a+b, then an error is indicated ("error[1]" is set to    */
/*  7).  If a+b is odd, p divides a+b, and p**2 divides a+b, then an error   */
/*  is indicated ("error[1]" is set to 7).                                   */
/*									     */
/*  If a is even, p does not divide a, and (a/2)**(p-1) is not congruent to  */
/*  1 modulo p**3, then an error is indicated ("error[3]" is set to 9).      */
/*  b is treated similarly.						     */
/*									     */
/*  If a is odd, p divides a+b, and a**(p-1) is not congruent to 1 modulo    */
/*  p**3, then an error is indicated ("error[3]" is set to 9).  b and a-b    */
/*  are treated similarly.						     */
/*									     */
/*  If p**3 divides a, b, a-b or p**4 divides a+b, 2 does not divide a, and  */
/*  a**(p-1) is not congruent to 1 modulo p**3, then an error is indicated   */
/*  ("error[2]" is set to 8). b and a-b are treated similarly.  If p**3      */
/*  divides a, b, or a-b or p**4 divides a+b, 2 does not divide a+b, p does  */
/*  not divide a+b, and (a+b)**(p-1) is not congruent to 1 modulo p**3,      */
/*  then an error is indicated.  If p**3 divides a, b, or a-b or p**4	     */
/*  divides a+b, 2 does not divide a+b, p divides a+b, and [(a+b)/p]**(p-1)  */
/*  is not congruent to 1 modulo p**3, then an error is indicated.	     */
/*									     */
/*  Proposition (28) can be verified by examining the output and determining */
/*  if p**3 divides a', b', or a'-b' or p**4 divides a'+b' for some          */
/*  representation [((a')**p+(b')**p)/(a'+b')] of T**p.  Use "prop28c.c"     */
/*  to group the representations of T**p and factor a, b, a-b, and a+b (the  */
/*  output of this program is input to "prop28c.c").                         */
/*									     */
/*  Note:  If p**4 divides a+b and [(a**p+b**p)/(a+b)] has only one distinct */
/*  factor, 2 does not divide a already implies that a**(p-1) is congruent   */
/*  to 1 modulo p**3 (similarly for b and a-b).  If p**3 divides a, b, or    */
/*  a-b and [(a**p+b**p)/(a+b)] has only one distinct prime factor, then     */
/*  p**3 divides a, b, or a-b for the other two representations (although    */
/*  this is probably not true for p>3).  The stipulation that p**3 divide    */
/*  a', b', or a'-b' or p**4 divide a'+b' is then really only useful when    */
/*  p**4 divides a'+b' (which is enough).  If p*T**p has two distinct prime  */
/*  factors, then [(a+b)/p]**(p-1) is not congruent to 1 modulo p**3 when    */
/*  p**4 divides a'+b' for the other group of three representations.         */
/*									     */
/*  Note:  If 2 does not divide a, the condition that p**3 divide a, b, or   */
/*  a-b or p**4 divide a+b is needed in case 2*p divides b or a-b	     */
/*  (similarly for b and a-b).						     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0b.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
void dummy(unsigned int a, unsigned int b, unsigned int c);
void sum(unsigned int *addend, unsigned int *augend);
void differ(unsigned int *minuend, unsigned int *subtrahend);
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void quotient(unsigned int *a, unsigned int *b, unsigned int);
int main ()
{
unsigned int p=3;	   // input prime
unsigned int dbeg=10000;   // starting "a" value
unsigned int dend=1;	   // ending "a" value
//unsigned int stop=0x2e831;
unsigned int sumdif=1;	   // select [(a**p+b**p)/(a+b)] if "sumdif" is non-
			   // zero, or [(a**p-b**p)/(a-b)] otherwise

extern unsigned int table3[];
extern unsigned int output[];
extern unsigned int error[];
unsigned int t3size=848;
unsigned int outsiz=1999;
unsigned int n=0;
unsigned int d,e,a,b,temp;
unsigned int i,j,k,l,lp,m,wrap;
unsigned int ps,pc,pf,qflag,flag,split;
unsigned int S[2],T[2],V[2],X[3];
double croot2,croot4,halfcr4;
FILE *Outfp;
Outfp = fopen("out28a.dat","w");
croot2=1.259921;
croot4=1.587401;
halfcr4=croot4*((double)(0.5));
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p;
pc=ps*p;
pf=pc*p;
error[0]=0;	// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
wrap=1;
for (d=dbeg; d>=dend; d--) {
   if (wrap>99) {
      printf("a=%d \n",d);
      wrap=1;
      }
   else
      wrap=wrap+1;
   for (e=d-1; e>0; e--) {
      dummy(d,e,2);
//    if (e!=stop) continue;
/*******************************/
/*  check for common factors   */
/*******************************/
      if((d==(d/2)*2)&&(e==(e/2)*2)) continue;
      if((d==(d/3)*3)&&(e==(e/3)*3)) continue;
      if((d==(d/5)*5)&&(e==(e/5)*5)) continue;
      if((d==(d/7)*7)&&(e==(e/7)*7)) continue;
/***********************/
/*  Euclidean G.C.D.   */
/***********************/
      a=d;
      b=e;
      if (b>a) {
	 temp=a;
	 a=b;
	 b=temp;
	 }
loop: temp = a - (a/b)*b;
      a=b;
      b=temp;
      if (b!=0) goto loop;
      if (a!=1) continue;
/*******************************/
/* check if 2 and p are split  */
/********************************/
      split=1;
      if ((d/2)*2==d)
	 if ((d/p)*p==d)
	    split=0;
      if ((e/2)*2==e)
	 if ((e/p)*p==e)
	    split=0;
      if (((d-e)/2)*2==(d-e))
	 if (((d-e)/p)*p==(d-e))
	    split=0;
      if (((d+e)/2)*2==(d+e))
	 if (((d+e)/p)*p==(d+e))
	    split=0;
/**********************************************************/
/* check if p**3 divides a, b, or a-b or p**4 divides a+b */
/**********************************************************/
      qflag=0;
      if ((d/pc)*pc==d)
	qflag=1;
      if ((e/pc)*pc==e)
	qflag=1;
      if (sumdif==1) {
	if (((d+e)/pf)*pf==(d+e))
	  qflag=1;
	if (((d-e)/pc)*pc==(d-e))
	  qflag=1;
	}
      if (sumdif==0) {
	if (((d-e)/pf)*pf==(d-e))
	  qflag=1;
	if (((d+e)/pc)*pc==(d+e))
	  qflag=1;
	}
/************************************/
/*  compute (d**p + e**p)/(d + e)   */
/************************************/
      S[0]=0;
      S[1]=d;
      for (i=0; i<p-1; i++) {
	 bigprod(S[0], S[1], d, X);
	 S[0]=X[1];
	 S[1]=X[2];
	 }
      T[0]=0;
      T[1]=e;
      for (i=0; i<p-1; i++) {
	 bigprod(T[0], T[1], e, X);
	 T[0]=X[1];
	 T[1]=X[2];
	 }
      if (sumdif==1) {
	 sum(S, T);
	 temp=d+e;
	 if (((d+e)/p)*p==(d+e))
	    temp=temp*p;
	 quotient(T, S, temp);
	 }
      else {
	 differ(S, T);
	 temp=d-e;
	 if (((d-e)/p)*p==(d-e))
	    temp=temp*p;
	 quotient(T, S, temp);
	 }
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
      if (S[0]==0)
	 l = 32 - lmbd(1, S[1]);
      else
	 l = 64 - lmbd(1, S[0]);
      j=l-(l/3)*3;
      l=l/3;
      l = 1 << l;
      if (j==0)
	 lp=(int)(((double)(l))*halfcr4);
      if (j==1) {
	 lp=l;
	 l=(int)(((double)(l))*croot2);
	 }
      if (j==2){
	 lp=(int)(((double)(l))*croot2);
	 l=(int)(((double)(l))*croot4);
	 }
      lp=lp-1;
      l=l+1;
      if (l>table3[t3size-1]) {
	 error[0]=5;
	 goto bskip;
	 }
      else {
	 j=0;
	 for (i=0; i<t3size; i++) {
	    if (table3[i] < lp) j=i;
	    else break;
	    }
	 k=j;
	 for (i=j; i<t3size; i++) {
	    if (table3[i] < l) k=i;
	    else break;
	    }
	 }
      for (i=j; i<=k; i++) {
	 m=0;
	 l = table3[i];
	 quotient(S, V, l);
	 bigprod(V[0], V[1], l, X);
	 if ((S[0]!=X[1]) || (S[1]!=X[2])) continue;
aloop:	 S[0]=V[0];
	 S[1]=V[1];
	 m=m+1;
	 quotient(S, V, l);
	 bigprod(V[0], V[1], l, X);
	 if ((S[0]==X[1]) && (S[1]==X[2])) goto aloop;
	 if ((m/3)*3!=m)
	    goto askip;
	 else
	    break;
	 }
      if ((m/3)*3!=m) continue;
      if ((S[0]!=0) || (S[1]!=1)) continue;
//
// p**2
//
      if ((d/2)*2==d) {
	 if ((d/p)*p==d) {
	    if ((d/ps)*ps!=d) {
	       error[1]=7;
	       goto bskip;
	       }
	    }
	 }
      if ((e/2)*2==e) {
	 if ((e/p)*p==e) {
	    if ((e/ps)*ps!=e) {
	       error[1]=7;
	       goto bskip;
	       }
	    }
	 }
      if (sumdif!=0) {
	 if (((d-e)/2)*2==(d-e)) {
	    if (((d-e)/p)*p==(d-e)) {
	       if (((d-e)/ps)*ps!=(d-e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 }
      else {
	 if (((d+e)/2)*2==(d+e)) {
	    if (((d+e)/p)*p==(d+e)) {
	       if (((d+e)/ps)*ps!=(d+e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 }
      if (sumdif!=0) {
	 if (((d+e)/2)*2==(d+e)) {
	    if (((d+e)/p)*p==(d+e)) {
	       if (((d+e)/pc)*pc!=(d+e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 if (((d+e)/2)*2!=(d+e)) {
	    if (((d+e)/p)*p==(d+e)) {
	       if (((d+e)/ps)*ps==(d+e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 }
      else {
	 if (((d-e)/2)*2==(d-e)) {
	    if (((d-e)/p)*p==(d-e)) {
	       if (((d-e)/pc)*pc!=(d-e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 if (((d-e)/2)*2!=(d-e)) {
	    if (((d-e)/p)*p==(d-e)) {
	       if (((d-e)/ps)*ps==(d-e)) {
		  error[1]=7;
		  goto bskip;
		  }
	       }
	    }
	 }
//
// p**3
//
      if (qflag!=0) {
	 if ((d/2)*2!=d) {
	    flag=0;
	    if (((d-1)/pc)*pc==(d-1))
	       flag=1;
	    if (((d+1)/pc)*pc==(d+1))
	       flag=1;
	    if (flag==0)
	       error[2]=8;
	    }
	 if ((e/2)*2!=e) {
	    flag=0;
	    if (((e-1)/pc)*pc==(e-1))
	       flag=1;
	    if (((e+1)/pc)*pc==(e+1))
	       flag=1;
	    if (flag==0)
	       error[2]=8;
	    }
	 if (sumdif!=0) {
	    if (((d-e)/2)*2!=(d-e)) {
	       flag=0;
	       if ((((d-e)-1)/pc)*pc==((d-e)-1))
		  flag=1;
	       if ((((d-e)+1)/pc)*pc==((d-e)+1))
		  flag=1;
	       if (flag==0)
		  error[2]=8;
	       }
	    if (((d+e)/2)*2!=(d+e)) {
	       if (((d+e)/p)*p==(d+e)) {
		  flag=0;
		  if (((((d+e)/p)-1)/pc)*pc==(((d+e)/p)-1))
		     flag=1;
		  if (((((d+e)/p)+1)/pc)*pc==(((d+e)/p)+1))
		     flag=1;
		  if (flag==0)
		     error[2]=8;
		  }
	       else {
		  flag=0;
		  if ((((d+e)-1)/pc)*pc==((d+e)-1))
		     flag=1;
		  if ((((d+e)+1)/pc)*pc==((d+e)+1))
		     flag=1;
		  if (flag==0)
		     error[2]=8;
		  }
	       }
	    }
	 else {
	    if (((d+e)/2)*2!=(d+e)) {
	       flag=0;
	       if ((((d+e)-1)/pc)*pc==((d+e)-1))
		  flag=1;
	       if ((((d+e)+1)/pc)*pc==((d+e)+1))
		  flag=1;
	       if (flag==0)
		  error[2]=8;
	       }
	    if (((d-e)/2)*2!=(d-e)) {
	       if (((d-e)/p)*p==(d-e)) {
		  flag=0;
		  if (((((d-e)/p)-1)/pc)*pc==(((d-e)/p)-1))
		     flag=1;
		  if (((((d-e)/p)+1)/pc)*pc==(((d-e)/p)+1))
		     flag=1;
		  if (flag==0)
		     error[2]=8;
		  }
	       else {
		  flag=0;
		  if ((((d-e)-1)/pc)*pc==((d-e)-1))
		     flag=1;
		  if ((((d-e)+1)/pc)*pc==((d-e)+1))
		     flag=1;
		  if (flag==0)
		     error[2]=8;
		  }
	       }
	    }
	 }
//
// "split" tests
//
      if ((d/2)*2==d) {
	 if ((d/p)*p!=d) {
	    flag=0;
	    if ((((d/2)-1)/pc)*pc==((d/2)-1))
	       flag=1;
	    if ((((d/2)+1)/pc)*pc==((d/2)+1))
	       flag=1;
	    if (flag==0)
	       error[3]=9;
	    }
	 }
      else {
	 if (sumdif==1) {
	    if (((d+e)/p)*p==(d+e)) {
	       flag=0;
	       if (((d-1)/pc)*pc==(d-1))
		  flag=1;
	       if (((d+1)/pc)*pc==(d+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 else {
	    if (((d-e)/p)*p==(d-e)) {
	       flag=0;
	       if (((d-1)/pc)*pc==(d-1))
		  flag=1;
	       if (((d+1)/pc)*pc==(d+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 }
      if ((e/2)*2==e) {
	 if ((e/p)*p!=e) {
	    flag=0;
	    if ((((e/2)-1)/pc)*pc==((e/2)-1))
	       flag=1;
	    if ((((e/2)+1)/pc)*pc==((e/2)+1))
	       flag=1;
	    if (flag==0)
	       error[3]=9;
	    }
	 }
      else {
	 if (sumdif==1) {
	    if (((d+e)/p)*p==(d+e)) {
	       flag=0;
	       if (((e-1)/pc)*pc==(e-1))
		  flag=1;
	       if (((e+1)/pc)*pc==(e+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 else {
	    if (((d-e)/p)*p==(d-e)) {
	       flag=0;
	       if (((e-1)/pc)*pc==(e-1))
		  flag=1;
	       if (((e+1)/pc)*pc==(e+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 }
      if (sumdif==1) {
	 if (((d-e)/2)*2!=(d-e)) {
	    if (((d+e)/p)*p==(d+e)) {
	       flag=0;
	       if ((((d-e)-1)/pc)*pc==((d-e)-1))
		  flag=1;
	       if ((((d-e)+1)/pc)*pc==((d-e)+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 }
      else {
	 if (((d+e)/2)*2!=(d+e)) {
	    if (((d-e)/p)*p==(d-e)) {
	       flag=0;
	       if ((((d+e)-1)/pc)*pc==((d+e)-1))
		  flag=1;
	       if ((((d+e)+1)/pc)*pc==((d+e)+1))
		  flag=1;
	       if (flag==0)
		  error[3]=9;
	       }
	    }
	 }
      if (n+2>outsiz) {
	 error[0]=6;
	 goto bskip;
	 }
      output[n]=d;
      output[n+1]=e;
      output[n+2]=l;
      output[n+3]=split;
      n=n+4;
askip:temp=0;
      }
   }
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d error2=%d error3=%d \n",error[0],
	error[1],error[2],error[3]);
fprintf(Outfp," count=%d \n",(n+1)/4);
for (i=0; i<(n+1)/4; i++)
   fprintf(Outfp," %#10x, %#10x, %#10x, %d, \n",output[4*i],output[4*i+1],
	   output[4*i+2],output[4*i+3]);
fclose(Outfp);
if ((error[1]!=0)||(error[2]!=0)||(error[3]!=0))
   printf(" error \n");
return(0);
}