/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/27/06 (dkc)							     */
/*									     */
/*  This C program determines if a, p*b, and a+b are pth powers w.r.t.	     */
/*  (a**p+b**p)/(a+b).	p**2 must divide a, b, a-b, or a+b.  Whether b, a-b, */
/*  and a+b are pth powers modulo p**2 when p**2 divides a is determined.    */
/*									     */
/*  Note:  The least-residues modulo p**2 table ("residue") is dependent on  */
/*  p.	Modify the look-up table accordingly.				     */
/*									     */
/*  If b, a-b, or a+b is not a pth power modulo p**2 when p**2 divides a,    */
/*  then an error is indicated ("error[1]" is set to a non-zero value).  If  */
/*  p and b are pth powers w.r.t. (a**p+b**p)/(a+b) and a-b is not a pth     */
/*  power w.r.t. (a**p+b**p)/(a+b), then an error is indicated.  If a-b is a */
/*  pth power w.r.t. (a**p+b**p)/(a+b) and p or b is not a pth power w.r.t.  */
/*  (a**p+b**p)/(a+b), then an error is indicated.			     */
/*									     */
/*  Note: Prime [(a^p+b^p)/(a+b)] (or [(a^p+b^p)/(a+b)] with small prime     */
/*  factors) are usually not checked.					     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table12.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void bigbigd(unsigned int *a, unsigned int *b);
void differ(unsigned int *a, unsigned int *b);
void dummy(unsigned int a, unsigned int b, unsigned int c);
void bigbigs(unsigned int *addend, unsigned int *augend);
void bigbigq(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f, unsigned int g);
void bigresx(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f);
void quotient(unsigned int *a, unsigned int *b, unsigned int c);
void newprod(unsigned int a, unsigned int b, unsigned int c, unsigned *d,
	     unsigned int e);
int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 2^32.
//	 The maximum "dbeg" value for p=5 is about 2^20.
//	 The maximum "dbeg" value for p=7 is about 2^14.
//	 The maximum "dbeg" value for p=11 is about 2^9.
//
unsigned int p=3;	  // input prime
unsigned int dbeg=180000;   // starting "a" value
unsigned int dend=1;	  // ending "a" value
//unsigned int stop=803;
unsigned int sumdif=1;	  // select [(a**p+b**p)/(a+b)] if "sumdif" is non-zero
			  // or [(a**p-b**p)/(a-b)] otherwise
unsigned int correct=0;
		// There is a small probability that (a**p+b**p)/(a+b) is not
		// completely factored if "correct" is not set to 1.
unsigned int residue[2]={1,8};			 // p=3
//unsigned int residue[4]={1,7,18,24};			   // p=5
//unsigned int residue[6]={1,18,19,30,31,48};		     // p=7
//unsigned int residue[10]={1,3,9,27,40,81,94,112,118,120};  // p=11

extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int tmpsav;
extern unsigned int count;
unsigned int pcount;
unsigned int qcount;
unsigned int maxsiz=100000;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int outsiz=499;
unsigned int save[16];	 // solutions array
unsigned int savsiz=15;  // size of solutions array minus one
unsigned int d,e,a,b,temp;
unsigned int i,j,k,l,m;
unsigned int flag,tflag,aflag,bflag,limit;
unsigned int S[4],T[4],U[4],V[4],W[4],X[4],Y[4],Z[4];
unsigned short rem;
unsigned int ps,iters,wrap;
unsigned int n=0;
double sqrt2=1.4142135;
unsigned int histom[16],histop[16];
FILE *Outfp;
Outfp = fopen("out38b.dat","w");
iters=0;
for (i=0; i<16; i++) {
   histom[i]=0;
   histop[i]=0;
   }
/*********************************/
/*  extend prime look-up table	 */
/*********************************/
tmpsiz=0;
for (i=0; i<tsize; i++) {
   j = (int)(table[i]);
   if (((j-1)/p)*p==(j-1)) {
      tmptab[tmpsiz] = j;
      tmpsiz=tmpsiz+1;
      }
   }
for (d=9975; d<(9971*9971); d++) {
   if (((d-1)/p)*p!=(d-1))
      continue;
   if(d==(d/2)*2) continue;
   if(d==(d/3)*3) continue;
   if(d==(d/5)*5) continue;
   if(d==(d/7)*7) continue;
   if(d==(d/11)*11) continue;
   if(d==(d/13)*13) continue;
   if(d==(d/17)*17) continue;
   if(d==(d/19)*19) continue;
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
   l = (int)(2.0 + sqrt((double)d));
   k=0;
   if (l>table[tsize-1]) {
      error[0]=1;
      goto bskip;
      }
   else {
      for (i=0; i<tsize; i++) {
	 if (table[i] < l) k=i;
	 else break;
	 }
      }
   flag=1;
   l=k;
   for (i=0; i<=l; i++) {
      k = table[i];
      if ((d/k)*k == d) {
	 flag=0;
	 break;
	 }
      }
   if (flag==1) {
      tmptab[tmpsiz]=d;
      tmpsiz = tmpsiz + 1;
      if (tmpsiz>=maxsiz)
	 break;
      }
   }
printf("size=%d, prime=%d \n",tmpsiz,tmptab[tmpsiz-1]);
tmpsav=tmpsiz;
limit=(tmptab[tmpsiz-1])>>16;
limit=limit*limit;
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p;
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
count=0;
pcount=0;
qcount=0;
wrap=0;
for (d=dbeg; d>=dend; d--) {
   if (p>5) {
      if (wrap>8) {
	 printf("d=%d \n",d);
	 wrap=0;
	 }
      else
	 wrap=wrap+1;
      }
   for (e=d-1; e>0; e--) {
//    if (e!=stop) continue;
/******************************************/
/* check if p^2 divides a, b, a+b, or a-b */
/******************************************/
      if ((d/ps)*ps==d)
	 goto dskip;
      if ((e/ps)*ps==e)
	 goto dskip;
      if (((d-e)/ps)*ps==(d-e))
	 goto dskip;
      if (((d+e)/ps)*ps!=(d+e))
	 continue;
/******************************************/
/*  check for common factors of d and e   */
/******************************************/
dskip: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 for Furtwangler conditions   */
/***************************************/
      tflag=1;
      if ((d/p)*p!=d) {
	 flag=0;
	 rem=d-(d/ps)*ps;
	 for (l=0; l<p-1; l++) {
	    if (rem==residue[l])
	       flag=1;
	    }
	 if (flag==0)
	    tflag=0;
	 }
      if ((e/p)*p!=e) {
	 flag=0;
	 rem=e-(e/ps)*ps;
	 for (l=0; l<p-1; l++) {
	    if (rem==residue[l])
	       flag=1;
	    }
	 if (flag==0)
	    tflag=0;
	 }
      if ((((d-e)/p)*p!=(d-e))&&((d+e)/p)*p!=(d+e)) {
	 flag=0;
	 rem=(d-e)-((d-e)/ps)*ps;
	 for (l=0; l<p-1; l++) {
	    if (rem==residue[l])
	       flag=1;
	    }
	 if (flag==0)
	    tflag=0;
	 flag=0;
	 rem=(d+e)-((d+e)/ps)*ps;
	 for (l=0; l<p-1; l++) {
	    if (rem==residue[l])
	       flag=1;
	    }
	 if (flag==0)
	    tflag=0;
	 }
/************************************/
/*  compute (d**p + e**p)/(d + e)   */
/************************************/
      for (i=0; i<16; i++)
	 save[i]=0;
      Y[0]=0;
      Y[1]=0;
      Y[2]=0;
      Y[3]=d;
      for (i=0; i<p-1; i++)
	 newprod(Y[1],Y[2],Y[3],Y,d);
      Z[0]=0;
      Z[1]=0;
      Z[2]=0;
      Z[3]=e;
      for (i=0; i<p-1; i++)
	 newprod(Z[1],Z[2],Z[3],Z,e);
      if (sumdif==1) {
	 bigbigs(Y, Z);
	 temp=d+e;
	 if (((d+e)/p)*p==(d+e))
	    temp=temp*p;
	 bigbigq(Z[0], Z[1], Z[2], Z[3], Y, 0, temp);
	 }
      else {
	 bigbigd(Y, Z);
	 temp=d-e;
	 if (((d-e)/p)*p==(d-e))
	    temp=temp*p;
	 bigbigq(Z[0], Z[1], Z[2], Z[3], Y, 0, temp);
	 }
      S[0]=Y[0];
      S[1]=Y[1];
      S[2]=Y[2];
      S[3]=Y[3];
      W[0]=S[0];
      W[1]=S[1];
      W[2]=S[2];
      W[3]=S[3];
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
      if (S[0]!=0) {
	 printf("S too big \n");
	 goto zskip;
	 }
      if (S[1]!=0)
	 l=96-lmbd(1, S[1]);
      else {
	 if (S[2]!=0)
	    l=64-lmbd(1, S[2]);
	 else
	    l=32-lmbd(1, S[3]);
	 }
      j=l-(l/2)*2;
      l=l/2;
      if (l>28) {
	 flag=1;
	 k=tmpsiz-1;
	 goto ajump;
	 }
      l=1<<l;
      if (j==1)
	 l=(int)(((double)(l))*sqrt2);
      l=l+1;
      flag=0;
      if (l>tmptab[tmpsiz-1]) {
	 flag=1;
	 k=tmpsiz-1;
	 }
      else {
	 k=0;
	 for (i=0; i<tmpsiz; i++) {
	    if (tmptab[i] < l) k=i;
	    else break;
	    }
	 }
ajump:m=0;
      for (i=0; i<=k; i++) {
	 l = tmptab[i];
	 bigbigq(S[0],S[1],S[2],S[3],T,0,l);
	 V[0]=T[0];
	 V[1]=T[1];
	 V[2]=T[2];
	 V[3]=T[3];
	 if (T[0]!=0) {
	    printf("multiplier too large \n");
	    goto zskip;
	    }
	 newprod(T[1],T[2],T[3],X,l);
	 if ((S[0]!=X[0])||(S[1]!=X[1])||(S[2]!=X[2])||(S[3]!=X[3]))
	    continue;
aloop:	 S[0]=V[0];
	 S[1]=V[1];
	 S[2]=V[2];
	 S[3]=V[3];
	 save[m]=l;
	 if (m < savsiz) m=m+1;
	 else {
	    error[0]=3;
	    goto bskip;
	    }
	 bigbigq(S[0],S[1],S[2],S[3],T,0,l);
	 V[0]=T[0];
	 V[1]=T[1];
	 V[2]=T[2];
	 V[3]=T[3];
	 newprod(T[1],T[2],T[3],X,l);
	 if ((S[0]==X[0])&&(S[1]==X[1])&&(S[2]=X[2])&&(S[3]==X[3]))
	    goto aloop;
	 }
/***********************************************/
/*  output prime factors satisfying criterion  */
/***********************************************/
      if ((S[0]!=0)||(S[1]!=0)||((S[2]&0xc0000000)!=0))
	 goto askip;
      if (m==0)
	 goto askip;
      S[0]=S[2];
      S[1]=S[3];
      if ((S[0]!=0) || (S[1]!=1)) {
	 if ((flag==1) && (correct==1)) {
	    if (S[0]==0)
	       j = (32 - lmbd(1, S[1]));
	    else
	       j = (64 - lmbd(1, S[0]));
	    k=j-(j/2)*2;
	    j=j/2;
	    j = 1 << j;
	    if (k==1)
	       j=(int)(((double)(j))*sqrt2);
	    for (i=tmptab[tmpsiz-1]; i<j; i+=2*p) {
	       quotient(S, T, i);
	       bigprod(T[0], T[1], i, X);
	       if ((X[1]==S[0]) && (X[2]==S[1])) {
		  bigresx(0, (i-1)/p, 0, i, U, d);
		  if ((U[0]!=0)||(U[1]!=1))
		     goto askip;
		  if (T[0]<=limit) {   // largest prime in table is 0x126f5f
		     S[0]=T[0];      // for p=7
		     S[1]=T[1];
		     save[m]=i;
		     if (m < savsiz) m=m+1;
		     else {
			error[0]=3;
			goto bskip;
			}
		     goto cskip;
		     }
		  else {
		     error[0]=4;
		     goto bskip;
		     }
		  }
	       }
	    }
cskip:	 T[0]=0;
	 T[1]=1;
	 differ(S, T);
	 quotient(T, T, p);
	 aflag=1;
	 bflag=1;
	 bigresx(T[0], T[1], S[0], S[1], U, d);
	 if ((U[0]==0)&&(U[1]==1)) {
	    bigresx(T[0], T[1], S[0], S[1], U, p*e);
	    if ((U[0]!=0)||(U[1]!=1))
	       goto askip;
	    if (sumdif!=0)
	       bigresx(T[0], T[1], S[0], S[1], U, d+e);
	    else
	       bigresx(T[0], T[1], S[0], S[1], U, d-e);
	    if ((U[0]!=0)||(U[1]!=1))
	       goto askip;
	    if (sumdif!=0)
	       bigresx(T[0], T[1], S[0], S[1], U, d-e);
	    else
	       bigresx(T[0], T[1], S[0], S[1], U, d+e);
	    if ((U[0]!=0)||(U[1]!=1))
	       aflag=0;
	    bigresx(T[0], T[1], S[0], S[1], U, e);
	    if ((U[0]!=0)||(U[1]!=1))
	       bflag=0;
	    for (i=0; i<m; i++) {
	       l=save[i];
	       bigresx(0, (l-1)/p, 0, l, U, d);
	       if ((U[0]!=0)||(U[1]!=1))
		  goto askip;
	       bigresx(0, (l-1)/p, 0, l, U, p*e);
	       if ((U[0]!=0)||(U[1]!=1))
		  goto askip;
	       if (sumdif!=0)
		  bigresx(0, (l-1)/p, 0, l, U, d+e);
	       else
		  bigresx(0, (l-1)/p, 0, l, U, d-e);
	       if ((U[0]!=0)||(U[1]!=1))
		  goto askip;
	       if (sumdif!=0)
		  bigresx(0, (l-1)/p, 0, l, U, d-e);
	       else
		  bigresx(0, (l-1)/p, 0, l, U, d+e);
	       if ((U[0]!=0)||(U[1]!=1))
		  aflag=0;
	       bigresx(0, (l-1)/p, 0, l, U, e);
	       if ((U[0]!=0)||(U[1]!=1))
		  bflag=0;
	       }
	    if (aflag!=bflag) {
	       error[1]=11;
	       error[2]=d;
	       error[3]=e;
	       printf("aflag=%d, bflag=%d \n",aflag,bflag);
	       goto bskip;
	       }
	    histom[m+1]=histom[m+1]+1;
	    if ((m==1)&&(S[1]==save[0]))
	       histop[2]=histop[2]+1;
	    if ((m==2)&&(S[1]==save[0])&&(S[1]==save[1]))
	       histop[3]=histop[3]+1;
	    if ((m==3)&&(S[1]==save[0])&&(S[1]==save[1])&&(S[1]==save[2]))
	       histop[4]=histop[4]+1;
	    if (((m>0)&&(p!=3))||((m>1)&&(p==3))) {
	       printf("d=%d, e=%d, m=%d, flag=%d, f=%d, %d, %d, %d, %d, %d  %#010x %#010x \n",d,e,m+1,aflag,save[0],save[1],save[2],save[3],save[4],save[5],S[0],S[1]);
	       if (iters<5000) {
		  fprintf(Outfp,"d=%d, e=%d, m=%d, flag=%d, f=%d, %d, %d, %d, %d, %d  %#010x %#010x \n",d,e,m+1,aflag,save[0],save[1],save[2],save[3],save[4],save[5],S[0],S[1]);
		  iters=iters+1;
		  }
	       }
	    T[0]=0;
	    T[1]=0;
	    T[2]=S[0];
	    T[3]=S[1];
	    for (i=0; i<m; i++) {
	       l=save[i];
	       newprod(T[1],T[2],T[3],X,l);
	       T[0]=X[0];
	       T[1]=X[1];
	       T[2]=X[2];
	       T[3]=X[3];
	       }
	    if ((T[0]!=W[0])||(T[1]!=W[1])||(T[2]!=W[2])||(T[3]!=W[3])) {
	       error[0]=7;
	       printf("error: incorrect product \n");
	       goto bskip;
	       }
	    if (tflag==0) {
	       error[1]=7;
	       goto bskip;
	       }
	    }
	 else {
	    goto askip;
	    }
	 }
      else {
	 aflag=1;
	 bflag=1;
	 for (i=0; i<m; i++) {
	    l=save[i];
	    bigresx(0, (l-1)/p, 0, l, U, d);
	    if ((U[0]!=0)||(U[1]!=1))
	       goto askip;
	    bigresx(0, (l-1)/p, 0, l, U, p*e);
	    if ((U[0]!=0)||(U[1]!=1))
	       goto askip;
	    if (sumdif!=0)
	       bigresx(0, (l-1)/p, 0, l, U, d+e);
	    else
	       bigresx(0, (l-1)/p, 0, l, U, d-e);
	    if ((U[0]!=0)||(U[1]!=1))
	       goto askip;
	    if (sumdif!=0)
	       bigresx(0, (l-1)/p, 0, l, U, d-e);
	    else
	       bigresx(0, (l-1)/p, 0, l, U, d+e);
	    if ((U[0]!=0)||(U[1]!=1))
	       aflag=0;
	    bigresx(0, (l-1)/p, 0, l, U, e);
	    if ((U[0]!=0)||(U[1]!=1))
	       bflag=0;
	    }
	 if (aflag!=bflag) {
	    error[1]=11;
	    error[2]=d;
	    error[3]=e;
	    printf("aflag=%d, bflag=%d \n",aflag,bflag);
	    goto bskip;
	    }
	 histom[m]=histom[m]+1;
	 if ((m==2)&&(save[0]==save[1]))
	    histop[2]=histop[2]+1;
	 if ((m==3)&&(save[0]==save[1])&&(save[0]==save[2]))
	    histop[3]=histop[3]+1;
	 if ((m==4)&&(save[0]==save[1])&&(save[0]==save[2])&&(save[0]==save[3]))
	    histop[4]=histop[4]+1;
	 if (((m>1)&&(p!=3))||((m>2)&&(p==3))) {
	    printf("d=%d, e=%d, m=%d, flag=%d, f=%d, %d, %d, %d, %d, %d \n",d,e,m,aflag,save[0],save[1],save[2],save[3],save[4],save[5]);
	    if (iters<2000) {
	       fprintf(Outfp,"d=%d, e=%d, m=%d, flag=%d, f=%d, %d, %d, %d, %d, %d \n",d,e,m,aflag,save[0],save[1],save[2],save[3],save[4],save[5]);
	       iters=iters+1;
	       }
	    }
	 S[0]=0;
	 S[1]=0;
	 S[2]=0;
	 S[3]=1;
	 for (i=0; i<m; i++) {
	    l=save[i];
	    newprod(S[1],S[2],S[3],X,l);
	    S[0]=X[0];
	    S[1]=X[1];
	    S[2]=X[2];
	    S[3]=X[3];
	    }
	 if ((S[0]!=W[0])||(S[1]!=W[1])||(S[2]!=W[2])||(S[3]!=W[3])) {
	    error[0]=7;
	    printf("error: incorrect product \n");
	    goto bskip;
	    }
	 if (tflag==0) {
	    error[1]=7;
	    goto bskip;
	    }
	 }
askip:dummy(d,e,6);
      }
   }
bskip:
printf("\n");
fprintf(Outfp,"\n");
fprintf(Outfp," error0=%d error1=%d count=%d asave=%d bsave=%d \n",error[0],
		error[1],error[2],error[3],error[4]);
if (error[1]!=0)
   printf(" error=%d \n",error[1]);
printf("\n");
fprintf(Outfp,"\n");
for (i=0; i<16; i++)
   fprintf(Outfp,"i=%d h=%d, %d, \n",i,histom[i],histop[i]);
zskip:
return(0);
}