/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b) (when [(a**p+b**p)/(a+b)] is a pth power)	     */
/*  11/03/06 (dkc)							     */
/*									     */
/*  This program factors (a**p+b**p)/(a+b).  Use "table3b" (same as "table3")*/
/*  for a<1000000.  Use "table4b" (same as "table4" except for the range of  */
/*  a) for 1000000<a<2000000.  Use "table6b" (same as "table6") for 2000000< */
/*  a<2500000.	(Modify "insize" accordingly.)                               */
/*									     */
/*  Note: [(a**p+b**p)/(a+b)] can have only three distinct prime factors.    */
/*									     */
/*  The output is "a, b, (factor0<<20)|(factor1<<10)|factor2, (power0<<20)|  */
/*  (power1<<10)|power2" where "factor0", "factor1", and "factor2" are the   */
/*  distinct prime factors of [(a**p+b**p)/(a+b)].  "power0", "power1", and  */
/*  "power2" are the powers of these factors.                                */
/*									     */
/*  If a is even, p does not divide a, and (a/2)**(p-1) is not congruent to  */
/*  1 mod p**3, then an error is indicated ("error[1] 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 modulus   */
/*  p**3, then an error is indicated ("error[1]" 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 modulus p**3, then an error is indicated  */
/*  ("error[1]" 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 modulus 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 modulus p**3, then an error is indicated.	     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table6b.h"
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void quotient(unsigned int *a, unsigned int *b, unsigned int);
void bigbigs(unsigned int *a, unsigned int *b);
void bigbigd(unsigned int *a, unsigned int *b);
void bigbigq(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f, unsigned int g);

int main ()
{
unsigned int p=3;       // input prime

extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int input[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int tmpsav;
unsigned int c,ps,pc,pf;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int insiz=1068;  // table6b
//unsigned int insiz=2534;  // table4b
//unsigned int insiz=4284;    // table3b
unsigned int outsiz=5000*3;
unsigned int d,e,temp,qflag;
unsigned int i,j,k,l,m,aflag,bflag;
unsigned int flag,sflag,tflag,iters,sumdif;
int pflag;
unsigned int S[2],T[2],V[2],X[3],Y[4],Z[4];
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out28f.dat","w");
/*********************************/
/*  extend prime look-up table	 */
/*********************************/
for (i=0; i<tsize; i++) tmptab[i] = (int)(table[i]);
tmpsiz=tsize;
for (d=10001; d<160000; d++) {
   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)(100.0 + sqrt((double)d));
      k=0;
      if (l>table[tsize-1]) {
	 error[0]=1;
aspin:	 goto aspin;
	 }
      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;
	 }
   }
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
   ps=p*p;
   pc=ps*p;
   pf=pc*p;
   pflag=0;
   error[0]=0;	   // clear error array
   error[1]=0;
   error[2]=0;
   for (j=0; j<insiz; j++) {
zloop:if (pflag<2) {
         d=input[3*j];
         e=input[3*j+1];
         c=input[3*j+2];
         sumdif=1;
         }
      else {
         d=input[3*(j-1)+1];
         e=input[3*j+1];
         c=input[3*j+2];
         sumdif=0;
         }
      if (c!=3)
         goto askip;
      if (sumdif==1) {
	 if (((d+e)/p)*p==(d+e))
	    qflag=1;
	 else
	    qflag=0;
	 }
      else {
	 if (((d-e)/p)*p==(d-e))
	    qflag=1;
	 else
	    qflag=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)   */
/************************************/
      Y[0] = 0;
      Y[1] = 0;
      Y[2] = 0;
      Y[3] = d;
      for (i=0; i<p-1; i++) {
	 bigprod(Y[2], Y[3], d, X);
	 Y[1]=X[0];
	 Y[2]=X[1];
	 Y[3]=X[2];
	 }
      Z[0] = 0;
      Z[1] = 0;
      Z[2] = 0;
      Z[3] = e;
      for (i=0; i<p-1; i++) {
	 bigprod(Z[2], Z[3], e, X);
	 Z[1]=X[0];
	 Z[2]=X[1];
	 Z[3]=X[2];
	 }
      if (sumdif!=0)
         bigbigs(Y, Z);
      else
         bigbigd(Y, Z);
      if (sumdif!=0) {
         temp=d+e;
         if (((d+e)/p)*p==(d+e))
            temp=temp*p;
         }
      else {
         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[2];
      S[1]=Y[3];
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
      iters=0;
      for (i=0; i<tmpsiz; i++) {
	 m=0;
	 k = tmptab[i];
	 quotient(S, T, k);
	 V[0]=T[0];
	 V[1]=T[1];
	 bigprod(T[0], T[1], k, X);
         if ((S[0]!=X[1]) || (S[1]!=X[2]))
            continue;
	 if (((k-1)/ps)*ps!=(k-1))
	    goto askip;
aloop:	 S[0]=V[0];
	 S[1]=V[1];
	 m=m+1;
	 quotient(S, T, k);
	 V[0]=T[0];
	 V[1]=T[1];
	 bigprod(T[0], T[1], k, X);
	 if ((S[0]==X[1]) && (S[1]==X[2])) goto aloop;
	 if ((m/3)*3!=m) {
	    error[0]=2;
	    goto bskip;
	    }
	 iters=iters+1;
	 if (iters==3)
	    break;
	 aflag=sflag;
	 sflag=k;
	 bflag=tflag;
	 tflag=m;
	 }
      if ((S[0]!=0)||(S[1]!=1)) {
	 error[0]=3;
	 goto bskip;
	 }
      sflag=(k<<20)|(sflag<<10)|aflag;
      tflag=(m<<20)|(tflag<<10)|bflag;
/***********************************************/
/*  output prime factors satisfying criterion  */
/***********************************************/
      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[1]=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[1]=8;
	    }
	 if (sumdif==1) {
	    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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=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[1]=9;
	       }
	    }
	 }
      if (n+4>outsiz) {
	 error[0]=6;
	 goto cskip;
	 }
      output[n]=d;
      output[n+1]=e;
      output[n+2]=sflag;
      output[n+3]=tflag;
      n=n+4;
askip:if (pflag==2)
         pflag=-1;
      pflag+=1;
      if (pflag==2)
         goto zloop;
      }
goto cskip;
bskip:
error[1]=d;
error[2]=e;
cskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d \n",error[0],error[1]);
fprintf(Outfp," count=%d \n",(n+1)/4);
for (i=0; i<(n+1)/4; i++)
   fprintf(Outfp," %#10x, %#10x, %#10x, %#10x, \n",output[4*i],output[4*i+1],
	   output[4*i+2],output[4*i+3]);
fclose(Outfp);
if (error[1]!=0)
   printf(" error \n");
return(0);
}