/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/08/06 (dkc)							     */
/*									     */
/*  This C program determines if q is a pth power with respect to	     */
/*  (a**p + b**p)/(a+b).  q must divide a, b, a+b or a-b.  If p>3, whether   */
/*  p divides a when q divides a and q is not a pth power modulo p**2	     */
/*  is determined (simiarly for b).  Whether p divides a-b or a+b when q     */
/*  divides a-b or a+b when q is not a pth power modulo p**2 is determined.  */
/*  q should be a prime.						     */
/*									     */
/*  The output is "(a<<16)|b".  If p does not divide a when q divides a,     */
/*  then an error is indicated ("error[1]" is set to a non-zero value).      */
/*  b is treated similarly.  If p does not divide a-b or a+b when q divides  */
/*  a-b or a+b, then an error is indicated.				     */
/*									     */
/*  Note: The user must ascertain that "base" is not a pth power modulo p^2. */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table11.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 hugeprod(unsigned int a, unsigned int b, unsigned int c, unsigned d,
	      unsigned int *e, unsigned int f);
void bigbigq(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f, unsigned int g);
void shift16(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e);
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);

int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 5000.
//	 The maximum "dbeg" value for p=5 is about 1000.
//	 The maximum "dbeg" value for p=7 is about 250.
//	 The maximum "dbeg" value for p=11 is about 50.
//
unsigned int p=7;	     // input prime
unsigned int dbeg=200;	     // starting "a" value
unsigned int dend=1;	     // ending "a" value
//unsigned int stop=0x3c1;
unsigned int base=29;	    // q value
	      // 1, 8 are pth powers modulus p**2 for p=3
	      // 1, 7, 18, 24 are pth powers modulus p**2 for p=5
	      // 1, 18, 19, 30, 31, 48 are pth powers modulus p**2 for p=7
	      // 1, 3, 9, 27, 40, 81, 94, 112, 118, 120 are pth powers
	      // modulus p**2 for p=11
unsigned int sumdif=1;	 // select [(a**p+b**p)/(a+b)] when "sumdif" is non-
			 // zero, or [(a**p-b**p)/(a-b)] otherwise
unsigned int correct=1;  // 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 out=1;    // results copied to the output array when set

extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int tmpsav;
extern unsigned int count;
extern unsigned int pcount;
unsigned int maxsiz=15600;
unsigned int tsize=303;
unsigned int tmpsiz;
unsigned int outsiz=1999;
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,limit,bsave;
unsigned int S[2],T[2],U[2],X[3],Y[4],Z[4];
unsigned int BT[4],BV[4],BW[4],BX[4];
unsigned int n=0;
double sqrt2=1.4142135;
FILE *Outfp;
Outfp = fopen("out23a.dat","w");
/*********************************/
/*  extend prime look-up table	 */
/*********************************/
error[0]=0;
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=2007; d<4000000; 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;
      }
   }
tmpsav=tmpsiz;
limit=(tmptab[tmpsiz-1])>>16;
limit=limit*limit;
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
error[6]=0;
count=0;
pcount=0;
for (d=dbeg; d>=dend; d--) {
   for (e=d-1; e>0; e--) {
//    if (e!=stop) continue;
/******************************************/
/*  check for common factors of d and e   */
/******************************************/
      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 q divides d, e, d+e or d-e   */
/******************************************/
      bsave=0;
      if (p!=3) {
	 if ((d/base)*base==d) {
	    bsave=d;
	    goto zskip;
	    }
	 if ((e/base)*base==e) {
	    bsave=e;
	    goto zskip;
	    }
	 }
      if (((d+e)/base)*base==(d+e))
	 goto zskip;
      if (((d-e)/base)*base!=(d-e))
	 continue;
/************************************/
/*  compute (d**p + e**p)/(d + e)   */
/************************************/
zskip:Y[0]=0;
      Y[1]=0;
      Y[2]=0;
      Y[3]=d;
      for (i=0; i<p-1; i++)
	 hugeprod(Y[0], 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++)
	 hugeprod(Z[0], 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);
	 }
      BW[0]=Y[0];
      BW[1]=Y[1];
      BW[2]=Y[2];
      BW[3]=Y[3];
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
      if ((Y[0]==0)&&(Y[1]==0)) {
	 if (Y[2]==0)
	    l = 32 - lmbd(1, Y[3]);
	 else
	    l = 64 - lmbd(1, Y[2]);
	 j=l-(l/2)*2;
	 l=l/2;
	 l = 1 << l;
	 if (j==1)
	    l=(int)(((double)(l))*sqrt2);
	 l=l+1;
	 }
      else
	 l=0x7fffffff;
      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;
	    }
	 }
      m=0;
      for (i=0; i<=k; i++) {
	 l = tmptab[i];
	 bigbigq(Y[0], Y[1], Y[2], Y[3], BT, 0, l);
	 BV[0]=BT[0];
	 BV[1]=BT[1];
	 BV[2]=BT[2];
	 BV[3]=BT[3];
	 if (l>65535) {
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BX, l>>16);
	    shift16(BX[0], BX[1], BX[2], BX[3], BX);
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l&65535);
	    bigbigs(BX, BT);
	    }
	 else
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l);
	 if ((Y[0]!=BT[0])||(Y[1]!=BT[1])||(Y[2]!=BT[2])||(Y[3]!=BT[3])) continue;
	 bigresx(0, (l-1)/p, 0, l, U, base);
	 if ((U[0]!=0)||(U[1]!=1)) {
	    goto askip;
	    }
aloop:	 Y[0]=BV[0];
	 Y[1]=BV[1];
	 Y[2]=BV[2];
	 Y[3]=BV[3];
	 save[m]=l;
	 if (m < savsiz) m=m+1;
	 else {
	    error[0]=3;
	    goto bskip;
	    }
	 bigbigq(Y[0], Y[1], Y[2], Y[3], BT, 0, l);
	 BV[0]=BT[0];
	 BV[1]=BT[1];
	 BV[2]=BT[2];
	 BV[3]=BT[3];
	 if (l>65535) {
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BX, l>>16);
	    shift16(BX[0], BX[1], BX[2], BX[3], BX);
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l&65535);
	    bigbigs(BX, BT);
	    }
	 else
	    hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l);
	 if ((Y[0]==BT[0])&&(Y[1]==BT[1])&&(Y[2]==BT[2])&&(Y[3]==BT[3])) goto aloop;
	 }
/***********************************************/
/*  output prime factors satisfying criterion  */
/***********************************************/
      if ((Y[0]!=0)||(Y[1]!=0))
	 continue;
      if (Y[2]>0x3fffffff)
	 continue;
      S[0]=Y[2];
      S[1]=Y[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, base);
		  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);
	 bigresx(T[0], T[1], S[0], S[1], U, base);
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (n>outsiz) {
	       error[0]=6;
	       goto bskip;
	       }
	    if (m>0)
	       printf("d=%d, e=%d, m=%d \n",d,e,m+1);
	    output[n]=((int)(d) << 16) | (int)(e);
	    BT[0]=0;
	    BT[1]=0;
	    BT[2]=S[0];
	    BT[3]=S[1];
	    for (i=0; i<m; i++) {
	       l=save[i];
	       if (l>65535) {
		  hugeprod(BT[0], BT[1], BT[2], BT[3], BX, l>>16);
		  shift16(BX[0], BX[1], BX[2], BX[3], BX);
		  hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l&65535);
		  bigbigs(BX, BT);
		  }
	       else
		  hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l);
	       }
	    if ((BT[0]!=BW[0])||(BT[1]!=BW[1])||(BT[2]!=BW[2])||(BT[3]!=BW[3])) {
	       error[0]=7;
	       goto bskip;
	       }
	    if (bsave!=0) {
	       if ((bsave/p)*p!=bsave) {
		  error[1]=12;
		  error[2]=d;
		  error[3]=e;
		  goto bskip;
		  }
	       }
	    else {
	       if ((((d+e)/p)*p!=(d+e))&&(((d-e)/p)*p!=(d-e))) {
		  error[1]=12;
		  error[2]=d;
		  error[3]=e;
		  goto bskip;
		  }
	       }
	    if (out!=0) {
	       n=n+1;
	       count=count+1;
	       }
	    else {
	       if (m>0) {
		  n=n+1;
		  count=count+1;
		  }
	       }
	    if (m>0) {
	       error[2]=error[2]+1;
	       if (out!=0)
		  m=error[2];
	       else
		  m=1;
	       if (m<100) {
		  error[2*m+3]=d;
		  error[2*m+4]=e;
		  }
	       }
	    }
	 else {
	    goto askip;
	    }
	 }
      else {
	 if (n>outsiz) {
	    error[0]=6;
	    goto bskip;
	    }
	 if (m>1)
	    printf("d=%d, e=%d, m=%d \n",d,e,m);
	 output[n]=((int)(d) << 16) | (int)(e);
	 BT[0]=0;
	 BT[1]=0;
	 BT[2]=0;
	 BT[3]=1;
	 for (i=0; i<m; i++) {
	    l=save[i];
	    if (l>65535) {
	       hugeprod(BT[0], BT[1], BT[2], BT[3], BX, l>>16);
	       shift16(BX[0], BX[1], BX[2], BX[3], BX);
	       hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l&65535);
	       bigbigs(BX, BT);
	       }
	    else
	       hugeprod(BT[0], BT[1], BT[2], BT[3], BT, l);
	    }
	 if ((BT[0]!=BW[0])||(BT[1]!=BW[1])||(BT[2]!=BW[2])||(BT[3]!=BW[3])) {
	    error[0]=7;
	    goto bskip;
	    }
	 if (bsave!=0) {
	    if ((bsave/p)*p!=bsave) {
	       error[1]=12;
	       error[2]=d;
	       error[3]=e;
	       goto bskip;
	       }
	    }
	 else {
	    if ((((d+e)/p)*p!=(d+e))&&(((d-e)/p)*p!=(d-e))) {
	       error[1]=12;
	       error[2]=d;
	       error[3]=e;
	       goto bskip;
	       }
	    }
	 if (out!=0) {
	    n=n+1;
	    count=count+1;
	    }
	 else {
	    if (m>1) {
	       n=n+1;
	       count=count+1;
	       }
	    }
	 if (m>1) {
	    error[2]=error[2]+1;
	    if (out!=0)
	       m=error[2];
	    else
	       m=1;
	    if (m<100) {
	       error[2*m+3]=d;
	       error[2*m+4]=e;
	       }
	    }
	 }
askip:dummy(d,e,6);
      }
   }
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d count=%d asave=%d bsave=%d \n",error[0],
		error[1],error[2],error[3],error[4]);
fprintf(Outfp," count=%d \n",n-1);
if (n!=0) {
   for (i=0; i<n-1; i++)
      fprintf(Outfp," %#10x \n",output[i]);
   }
fclose(Outfp);
if (error[1]!=0)
   printf(" error \n");
return(0);
}