/*****************************************************************************/
/*									     */
/*  FACTOR (a**p-b**p)/(a-b) (when [(a**p+b**p)/(a+b)] is a pth power)	     */
/*  02/25/14 (dkc)							     */
/*									     */
/*  Program determines if a, b, a-b, a+b, pa, pb, p(a-b), p(a+b), p**2*a,    */
/*  p**2*b, p**2*(a-b), p**2*(a+b), 2, p, 2p, or p/2 are pth powers modulo   */
/*  the prime factors of [(a**p-b**p)/(a-b)].  Use "table4a" for a<=2000000. */
/*  Use "table6a" for 2000000<a<=2500000.  Use "table7a" for a<2500000<=     */
/*  3000000.  Modify "insize" accordingly.                                   */
/*									     */
/*  The output is "a, b, (code<<16)|number of prime factors, code0||code1,   */
/*  code2||code3, code4||code5, code6, "and" of code0,code1,...,code6".      */
/*  "code" values are 0, 1, 2, or 3 if p divides a, b, a-b, or a+b           */
/*  respectively.  code0, code1,...,code6 are codes for the different prime  */
/*  factors of [(a**p-b**p)/(a-b)].  A corresponding bit is set in the code  */
/*  if p/2, 2*p, 2, p, a, b, a-b, a+b, pa, pb, p(a-b), p(a+b), p**2*a,	     */
/*  p**2*b, p**2*(a-b), or p**2*(a+b) is a pth power modulo a prime factor   */
/*  of [(a**p-b**p)/(a-b)].						     */
/*									     */
/*  If 2*p does not divide a, b, a-b, or a+b, and 2 is a pth power modulo    */
/*  a prime factor of [(a**p-b**p)/(a-b)] not of the form p**2*k+1, then     */
/*  an error is indicated ("error[2]" is set to 8).  (Set "split" to 1.)     */
/*  This is part of Proposition (4).					     */
/*									     */
/*  If 2*p divides a or b, [(a**p-b**p)/(a-b)] has two distinct prime	     */
/*  factors, neither of these distinct prime factors is of the form	     */
/*  p**2*k+1, and [(a**p-b**p)/(a-b)] is not of the form p**2*k+1, then      */
/*  an error is indicated ("error[2]" is set to 8).  This is part of         */
/*  Proposition (15).  The rest of Proposition (15) can be verified by	     */
/*  examining the codes generated when "numfac" equals 2, "split" equals 0,  */
/*  "psflag" equals 0, and "code" equals 0 or 1.                             */
/*									     */
/*  Proposition (8) can be verified by examining the codes generated when    */
/*  "psflag" equals 1.  (The codes should also be examined when "mixed"      */
/*  forms of prime factors of [(a**p-b**p)/(a-b)] are allowed.)  (Set	     */
/*  "split" to 2.)                                                           */
/*									     */
/*  Proposition (9) can be verified by examining the codes generated when    */
/*  "psflag" equals 0.  (The codes should also be examined when "mixed"      */
/*  forms of prime factors of [(a**p-b**p)/(a-b)] are allowed.)  (Set	     */
/*  "split" to 1.)                                                           */
/*									     */
/*  Proposition (3) can be verified by examining the "and" of the codes      */
/*  when "split" equals 1 and "code" equals 2 or 3.  The bits corresponding  */
/*  to 2 and p should be examined.					     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table7a.h"
void dummy(unsigned int a, unsigned int b, unsigned int c, unsigned int d);
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 bigresx(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f);
void bigbigs(unsigned int *a, unsigned int *b);
void bigbigd(unsigned int *a, unsigned int *b);
void bigbigq(unsigned int a0, unsigned int a1, unsigned int a2,
	     unsigned int a3, unsigned int *quotient, unsigned int d2,
	     unsigned int d3);
unsigned int lmbd(unsigned int mode, unsigned int a);
void differ(unsigned int *a, unsigned int *b);

int main ()
{
unsigned int p=3; // input prime
unsigned int numfac=3; // number of distinct prime factors of (a**p-b**p)/(a-b)
unsigned int split=2;  // if set to 0, don't allow 2 and p to "split"
		       // if set to 1, only allow "split" 2 and p
		       // otherwise, allow both
unsigned int psflag=2; // if set to 1, factors must be of the form p**2*k+1
		       // if set to 0, factors must not be of the form p**2*k+1
		       // otherwise, factors can be of mixed types
unsigned int offset=0;	   // offset into the input look-up table (must be
			   // even and less than "insize").  Used to reduce
			   // execution time.

extern unsigned int input[];
extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int count;
unsigned int insize=1038;  // table7a
//unsigned int insize=1068;   // table6a
//unsigned int insize=6818;   // table4a
unsigned int tsize=1228;
unsigned int maxsiz=200000;
unsigned int tmpsiz;
unsigned int outsiz=11000*7;
unsigned int save[20];  // solutions array
unsigned int savsiz=19;  // size of solutions array minus one
unsigned int d,e,c;
unsigned int h,i,j,k,l,m,iters;
unsigned int flag,temp,sumdif,ps,limit;
unsigned int R[2],S[2],T[2],U[2],V[2],W[2],X[3],Y[4],Z[4],Up[2];
int yflag,zflag;
unsigned int pflag,qflag,rflag,sflag,tflag,uflag,vflag,wflag,xflag,wrap;
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out8.dat","w");
if (numfac>7) {
   printf("number of factors too large (not enough memory allocated) \n");
   goto zskip;
   }
/*********************************/
/*  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=10001; d<10000000; 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)   */
/***********************************/
error[0]=0;	// clear error array
error[1]=0;
error[2]=0;
ps=p*p;
zflag=0;
count=0;
wrap=0;
for (h=offset; h<insize; h++) {
   if (wrap>8) {
      printf("count=%d \n",h+1);
      wrap=0;
      }
   else
      wrap=wrap+1;
zloop:
   if (zflag<2) {
      d=input[3*h];
      e=input[3*h+1];
      c=input[3*h+2];
      sumdif=0;
      }
   else {
      d=input[3*(h-1)+1];
      e=input[3*h+1];
      c=input[3*h+2];
      sumdif=1;
      }
//
// check for "split" 2 and p
//
      qflag=2;
      if ((d/2)*2==d) {
	 if (sumdif==0) {
	    if (((d+e)/p)*p==(d+e))
	       qflag=0;
	    }
	 else {
	    if (((d-e)/p)*p==(d-e))
	       qflag=0;
	    }
	 }
      if ((e/2)*2==e) {
	 if (sumdif==0) {
	    if (((d+e)/p)*p==(d+e))
	       qflag=1;
	    }
	 else {
	    if (((d-e)/p)*p==(d-e))
	       qflag=1;
	    }
	 }
      if ((split==0)&&(qflag!=2))
	 goto askip;
      if ((split==1)&&(qflag==2))
	 goto askip;
//
// check if p divides a, b, a-b, or a+b
//
      if ((d/p)*p==d)
	 yflag=0;
      if ((e/p)*p==e)
	 yflag=1;
      if (((d+e)/p)*p==(d+e)) {
	 if (sumdif==0)
	    yflag=3;
	 else
	    yflag=2;
	 }
      if (((d-e)/p)*p==(d-e)) {
	 if (sumdif==0)
	    yflag=2;
	 else
	    yflag=3;
	 }
/************************************/
/*  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==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[2];
      S[1]=Y[3];
      W[0]=S[0];
      W[1]=S[1];
/************************************************/
/*  look for prime factors using look-up table	*/
/************************************************/
      if (S[0]==0) {
	 l = (33 - lmbd(1, S[1]))/2;
	 l = 1 << l;
	 }
      else {
	 l = (65 - lmbd(1, S[0]))/2;
	 l = 1 << l;
	 }
      k=0;
      if (l>tmptab[tmpsiz-1]) {
	 flag=1;
	 k=tmpsiz-1;
	 }
      else {
	 flag=0;
	 for (i=0; i<tmpsiz; i++) {
	    if (tmptab[i] < l) k=i;
	    else break;
	    }
	 }
      l=k;
      iters=0;
      rflag=0xffffffff;
      sflag=0xffffffff;
      tflag=0xffffffff;
      uflag=0xffffffff;
      vflag=0xffffffff;
      wflag=0xffffffff;
      xflag=0xffffffff;
      m=0;
      for (i=0; i<=l; i++) {
	 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 (psflag==1) {
            if (((k-1)/ps)*ps!=(k-1))
               goto askip;
            }
	 if (psflag==0) {
	    if (((k-1)/ps)*ps==(k-1))
               goto askip;
            }
	 pflag=0;
	 if (((k-1)/ps)*ps==(k-1))
	    pflag=1;
	 if (iters<numfac) {
	    rflag=0;
	    bigresx(0, (k-1)/p, 0, k, U, d);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+2048;
	    bigresx(0, (k-1)/p, 0, k, U, e);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+1024;
	    bigresx(0, (k-1)/p, 0, k, U, d-e);
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+512;
	       else
		  rflag=rflag+256;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, d+e);
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+256;
	       else
		  rflag=rflag+512;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, p*d);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+128;
	    bigresx(0, (k-1)/p, 0, k, U, p*e);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+64;
	    bigresx(0, (k-1)/p, 0, k, U, p*(d-e));
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+32;
	       else
		  rflag=rflag+16;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, p*(d+e));
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+16;
	       else
		  rflag=rflag+32;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, p*p*d);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+8;
	    bigresx(0, (k-1)/p, 0, k, U, p*p*e);
	    if ((U[0]==0)&&(U[1]==1))
	       rflag=rflag+4;
	    bigresx(0, (k-1)/p, 0, k, U, p*p*(d-e));
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+2;
	       else
		  rflag=rflag+1;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, p*p*(d+e));
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (zflag!=2)
		  rflag=rflag+1;
	       else
		  rflag=rflag+2;
	       }
	    bigresx(0, (k-1)/p, 0, k, U, p);
	    Up[0]=U[0];
	    Up[1]=U[1];
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (pflag!=0) {
		  if (rflag!=0xfff)
		     error[2]=10;
		  }
	       rflag=rflag+4096;
	       }
	    else {
	       if (pflag!=0) {
		  if (yflag<3) {
		     if (rflag!=0xd20)
			error[2]=11;
		     }
		  else {
		     if (rflag!=0x0d2)
			error[2]=11;
		     }
		  }
	       else {
		  if (qflag==2) {
		     if (yflag==2) {
			if ((rflag!=0x1a4)&&(rflag!=0x168))
			   error[2]=12;
			}
		     if (yflag==3) {
			if ((rflag!=0x816)&&(rflag!=0x41a))
			   error[2]=12;
			}
		     }
		  }
	       }
	    R[0]=U[0];
	    R[1]=U[1];
	    bigresx(0, (k-1)/p, 0, k, U, 2);
	    if ((U[0]==Up[0])&&(U[1]==Up[1]))
	       rflag=rflag+32768;
	    if ((U[0]==0)&&(U[1]==1)) {
	       rflag=rflag+8192;
	       if (split==1) {
		  if (((k-1)/ps)*ps!=(k-1))
		     error[2]=8;
		  }
	       }
	    bigresx(0, (k-1)/p, 0, k, U, 2*p);
	    if ((U[0]==0)&&(U[1]==1)) {
	       if (pflag==0) {
		  if ((R[0]!=0)||(R[1]!=1)) {
		     if (qflag==2) {
			if (yflag==0) {
			   if (rflag!=0x816)
			      error[2]=15;
			   }
			if (yflag==1) {
			   if (rflag!=0x41a)
			      error[2]=15;
			   }
			}
		     else {
			if (qflag==0) {
			   if (rflag!=0x681)
			      error[2]=16;
			   }
			if (qflag==1) {
			   if (rflag!=0xa41)
			      error[2]=16;
			   }
			}
		     }
		  }
	       rflag=rflag+16384;
	       }
	    else {
	       if (pflag==0) {
		  if ((R[0]!=0)||(R[1]!=1)) {
		     if (qflag==2) {
			if (yflag==0) {
			   if (rflag!=0x8a41)
			      error[2]=13;
			   }
			if (yflag==1) {
			   if (rflag!=0x8681)
			      error[2]=13;
			   }
			}
		     else {
			if (qflag==0) {
			   if (rflag!=0x81a4)
			      error[2]=14;
			   }
			if (qflag==1) {
			   if (rflag!=0x8168)
			      error[2]=14;
			   }
			}
		     }
		  }
	       }
	    }
aloop:	 S[0]=V[0];
	 S[1]=V[1];
	 save[m]=k;
	 if (m < savsiz) m=m+1;
	 else {
	    error[0]=3;
	    goto bskip;
	    }
	 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;
	 iters=iters+1;
	 if (iters<numfac) {
	     xflag=wflag;
	     wflag=vflag;
	     vflag=uflag;
	     uflag=tflag;
	     tflag=sflag;
	     sflag=rflag;
	     }
	 else {
	    if ((S[0]!=0)||(S[1]!=1))
	       goto askip;
	    }
	 }
/***********************************************/
/*  output prime factors satisfying criterion  */
/***********************************************/
      if ((S[0]!=0) || (S[1]!=1)) {
	 if (flag==1) {
	    if (S[0]==0) {
	       j = (33 - lmbd(1, S[1]))/2;
	       j = 1 << j;
	       }
	    else {
	       j = (65 - lmbd(1, S[0]))/2;
	       j = 1 << j;
	       }
	    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])) {
		  if (psflag==1) {
                     if (((i-1)/ps)*ps!=(i-1))
                        goto askip;
                     }
		  if (psflag==0) {
		     if (((i-1)/ps)*ps==(i-1))
                        goto askip;
                     }
		  pflag=0;
		  if (((i-1)/ps)*ps==(i-1))
		     pflag=1;
		  if (iters<numfac) {
		     rflag=0;
		     bigresx(0, (i-1)/p, 0, i, U, d);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+2048;
		     bigresx(0, (i-1)/p, 0, i, U, e);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+1024;
		     bigresx(0, (i-1)/p, 0, i, U, d-e);
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+512;
			else
			   rflag=rflag+256;
			}
		     bigresx(0, (i-1)/p, 0, i, U, d+e);
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+256;
			else
			   rflag=rflag+512;
			}
		     bigresx(0, (i-1)/p, 0, i, U, p*d);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+128;
		     bigresx(0, (i-1)/p, 0, i, U, p*e);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+64;
		     bigresx(0, (i-1)/p, 0, i, U, p*(d-e));
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+32;
			else
			   rflag=rflag+16;
			}
		     bigresx(0, (i-1)/p, 0, i, U, p*(d+e));
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+16;
			 else
			   rflag=rflag+32;
			}
		     bigresx(0, (i-1)/p, 0, i, U, p*p*d);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+8;
		     bigresx(0, (i-1)/p, 0, i, U, p*p*e);
		     if ((U[0]==0)&&(U[1]==1))
			rflag=rflag+4;
		     bigresx(0, (i-1)/p, 0, i, U, p*p*(d-e));
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+2;
			else
			   rflag=rflag+1;
			}
		     bigresx(0, (i-1)/p, 0, i, U, p*p*(d+e));
		     if ((U[0]==0)&&(U[1]==1)) {
			if (zflag!=2)
			   rflag=rflag+1;
			else
			   rflag=rflag+2;
			}
		     bigresx(0, (i-1)/p, 0, i, U, p);
		     Up[0]=U[0];
		     Up[1]=U[1];
		     if ((U[0]==0)&&(U[1]==1)) {
			if (pflag!=0) {
			   if (rflag!=0xfff)
			      error[2]=10;
			   }
			rflag=rflag+4096;
			}
		     else {
			if (pflag!=0) {
			   if (yflag<3) {
			      if (rflag!=0xd20)
				 error[2]=11;
			      }
			   else {
			      if (rflag!=0x0d2)
				 error[2]=11;
			      }
			   }
			else {
			   if (qflag==2) {
			      if (yflag==2) {
				 if ((rflag!=0x1a4)&&(rflag!=0x168))
				    error[2]=12;
				 }
			      if (yflag==3) {
				 if ((rflag!=0x816)&&(rflag!=0x41a))
				    error[2]=12;
				 }
			      }
			   }
			}
		     R[0]=U[0];
		     R[1]=U[1];
		     bigresx(0, (i-1)/p, 0, i, U, 2);
		     if ((U[0]==Up[0])&&(U[1]==Up[1]))
			rflag=rflag+32768;
		     if ((U[0]==0)&&(U[1]==1)) {
			rflag=rflag+8192;
			if (split==1) {
			   if (((i-1)/ps)*ps!=(i-1))
			      error[2]=8;
			   }
			}
		     bigresx(0, (i-1)/p, 0, i, U, 2*p);
		     if ((U[0]==0)&&(U[1]==1)) {
			if (pflag==0) {
			   if ((R[0]!=0)||(R[1]!=1)) {
			      if (qflag==2) {
				 if (yflag==0) {
				    if (rflag!=0x816)
				       error[2]=15;
				    }
				 if (yflag==1) {
				    if (rflag!=0x41a)
				       error[2]=15;
				    }
				 }
			      else {
				 if (qflag==0) {
				    if (rflag!=0x681)
				       error[2]=16;
				    }
				 if (qflag==1) {
				    if (rflag!=0xa41)
				       error[2]=16;
				    }
				 }
			      }
			   }
			rflag=rflag+16384;
			}
		     else {
			if (pflag==0) {
			   if ((R[0]!=0)||(R[1]!=1)) {
			      if (qflag==2) {
				 if (yflag==0) {
				    if (rflag!=0x8a41)
				       error[2]=13;
				    }
				 if (yflag==1) {
				    if (rflag!=0x8681)
				       error[2]=13;
				    }
				 }
			      else {
				 if (qflag==0) {
				    if (rflag!=0x81a4)
				       error[2]=14;
				    }
				 if (qflag==1) {
				    if (rflag!=0x8168)
				       error[2]=14;
				    }
				 }
			      }
			   }
			}
		     }
		  iters=iters+1;
		  if (iters<numfac) {
		     xflag=wflag;
		     wflag=vflag;
		     vflag=uflag;
		     uflag=tflag;
		     tflag=sflag;
		     sflag=rflag;
		     }
		  if (T[0]<=limit) {	 // largest prime in table is 2751787
		     S[0]=T[0];
		     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:	 if ((S[0]==0)&&(S[1]==1)) {
	    if (iters!=numfac)
	       goto askip;
	    else
	       goto dskip;
	    }
	 if ((S[0]==0)&&(S[1]==save[m])) {
	    if (iters==numfac)
	       goto dskip;
	    else
	       goto askip;
	    }
	 if (iters!=(numfac-1))
	    goto askip;
	 if (psflag==1) {
            T[0]=0;
            T[1]=1;
            differ(S,T);
            quotient(T,T,ps);
            bigprod(T[0],T[1],ps,X);
            T[0]=0;
            T[1]=1;
            differ(S,T);
            if ((X[1]!=T[0])||(X[2]!=T[1]))
               goto askip;
            }
	 if (psflag==0) {
            T[0]=0;
            T[1]=1;
            differ(S,T);
            quotient(T,T,ps);
            bigprod(T[0],T[1],ps,X);
            T[0]=0;
            T[1]=1;
            differ(S,T);
	    if ((X[1]==T[0])&&(X[2]==T[1]))
               goto askip;
            }
	 pflag=0;
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 quotient(T,T,ps);
	 bigprod(T[0],T[1],ps,X);
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 if ((X[1]==T[0])&&(X[2]==T[1]))
	    pflag=1;
         T[0]=0;
	 T[1]=1;
	 differ(S, T);
	 quotient(T, T, p);
	 rflag=0;
	 bigresx(T[0], T[1], S[0], S[1], U, d);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+2048;
	 bigresx(T[0], T[1], S[0], S[1], U, e);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+1024;
	 bigresx(T[0], T[1], S[0], S[1], U, d-e);
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+512;
	    else
	       rflag=rflag+256;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, d+e);
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+256;
	    else
	       rflag=rflag+512;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, p*d);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+128;
	 bigresx(T[0], T[1], S[0], S[1], U, p*e);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+64;
	 bigresx(T[0], T[1], S[0], S[1], U, p*(d-e));
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+32;
	    else
	       rflag=rflag+16;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, p*(d+e));
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+16;
	    else
	       rflag=rflag+32;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, p*p*d);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+8;
	 bigresx(T[0], T[1], S[0], S[1], U, p*p*e);
	 if ((U[0]==0)&&(U[1]==1))
	    rflag=rflag+4;
	 bigresx(T[0], T[1], S[0], S[1], U, p*p*(d-e));
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+2;
	    else
	       rflag=rflag+1;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, p*p*(d+e));
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (zflag!=2)
	       rflag=rflag+1;
	    else
	       rflag=rflag+2;
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, p);
	 Up[0]=U[0];
	 Up[1]=U[1];
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (pflag!=0) {
	       if (rflag!=0xfff)
		  error[2]=10;
	       }
	    rflag=rflag+4096;
	    }
	 else {
	    if (pflag!=0) {
	       if (yflag<3) {
		  if (rflag!=0xd20)
		     error[2]=11;
		  }
	       else {
		  if (rflag!=0x0d2)
		     error[2]=11;
		  }
	       }
	    else {
	       if (qflag==2) {
		  if (yflag==2) {
		     if ((rflag!=0x1a4)&&(rflag!=0x168))
			error[2]=12;
		     }
		  if (yflag==3) {
		     if ((rflag!=0x816)&&(rflag!=0x41a))
			error[2]=12;
		     }
		  }
	       }
	    }
	 R[0]=U[0];
	 R[1]=U[1];
	 bigresx(T[0], T[1], S[0], S[1], U, 2);
	 if ((U[0]==Up[0])&&(U[1]==Up[1]))
	    rflag=rflag+32768;
	 if ((U[0]==0)&&(U[1]==1)) {
	    rflag=rflag+8192;
	    if (split==1) {
	       V[0]=0;
	       V[1]=1;
	       differ(S,V);
	       quotient(V,V,ps);
	       bigprod(V[0],V[1],ps,X);
	       V[0]=0;
	       V[1]=1;
	       differ(S,V);
	       if ((X[1]!=V[0])||(X[2]!=V[1]))
		  error[2]=8;
	       }
	    }
	 bigresx(T[0], T[1], S[0], S[1], U, 2*p);
	 if ((U[0]==0)&&(U[1]==1)) {
	    if (pflag==0) {
	       if ((R[0]!=0)||(R[1]!=1)) {
		  if (qflag==2) {
		     if (yflag==0) {
			if (rflag!=0x816)
			   error[2]=15;
			}
		     if (yflag==1) {
			if (rflag!=0x41a)
			   error[2]=15;
			}
		     }
		  else {
		     if (qflag==0) {
			if (rflag!=0x681)
			   error[2]=16;
			}
		     if (qflag==1) {
			if (rflag!=0xa41)
			   error[2]=16;
			}
		     }
		  }
	       }
	    rflag=rflag+16384;
	    }
	 else {
	    if (pflag==0) {
	       if ((R[0]!=0)||(R[1]!=1)) {
		  if (qflag==2) {
		     if (yflag==0) {
			if (rflag!=0x8a41)
			   error[2]=13;
			}
		     if (yflag==1) {
			if (rflag!=0x8681)
			   error[2]=13;
			}
		     }
		  else {
		     if (qflag==0) {
			if (rflag!=0x81a4)
			   error[2]=14;
			}
		     if (qflag==1) {
			if (rflag!=0x8168)
			   error[2]=14;
			}
		     }
		  }
	       }
	    }
dskip:	 if (n+6>outsiz) {
	    error[0]=6;
	    goto bskip;
	    }
	 output[n]=d;
	 output[n+1]=e;
	 output[n+2]=(yflag<<16)|(m+1);
	 output[n+3]=(xflag<<16)|wflag;
	 output[n+4]=(vflag<<16)|uflag;
	 output[n+5]=(tflag<<16)|sflag;
	 xflag=rflag&sflag&tflag&uflag&vflag&wflag&wflag;
	 output[n+6]=(rflag<<16)|xflag;
	 for (i=0; i<m; i++) {
	    bigprod(S[0], S[1], save[i], X);
	    S[0] = X[1];
	    S[1] = X[2];
	    }
	 if ((S[0]!=W[0]) || (S[1]!=W[1])) {
	    error[0]=7;
	    goto bskip;
	    }
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 quotient(T,T,ps);
	 bigprod(T[0],T[1],ps,X);
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 if ((X[1]!=T[0])||(X[2]!=T[1]))
	    error[1]+=1;
	 n=n+7;
	 count=count+1;
	 if ((numfac==2)&&(split==0)&&(psflag==0)) {
	    if ((yflag==0)||(yflag==1)) {
	       if (rflag!=sflag)
		  error[2]=8;
	       if ((rflag&0xf000)==0x2000)
		  error[2]=8;
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       quotient(T,T,ps);
	       bigprod(T[0],T[1],ps,X);
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       if ((X[1]!=T[0])||(X[2]!=T[1]))
		  error[2]=8;
	       }
	    else {
	       if ((rflag!=0xf333)&&(sflag!=0xf333)) {
		  if ((rflag&0xf000)!=0x2000)
		     error[2]=9;
		  if ((sflag&0xf000)!=0x2000)
		     error[2]=9;
		  }
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       quotient(T,T,ps);
	       bigprod(T[0],T[1],ps,X);
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       if ((X[1]!=T[0])||(X[2]!=T[1]))
		  error[2]=9;
	       }
	    }
	 }
//
// S[0]=0, S[1]=1
//
      else {
	 if (iters!=numfac)
	   goto askip;
	 if (n+6>outsiz) {
	    error[0]=6;
	    goto bskip;
	    }
	 output[n]=d;
	 output[n+1]=e;
	 output[n+2]=(yflag<<16)|m;
	 output[n+3]=(xflag<<16)|wflag;
	 output[n+4]=(vflag<<16)|uflag;
	 output[n+5]=(tflag<<16)|sflag;
	 xflag=rflag&sflag&tflag&uflag&vflag&wflag&wflag;
	 output[n+6]=(rflag<<16)|xflag;
         S[0]=0;
	 S[1]=1;
	 for (i=0; i<m; i++) {
	    bigprod(S[0], S[1], save[i], X);
	    S[0] = X[1];
	    S[1] = X[2];
	    }
	 if ((S[0]!=W[0]) || (S[1]!=W[1])) {
	    error[0]=7;
	    goto bskip;
	    }
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 quotient(T,T,ps);
	 bigprod(T[0],T[1],ps,X);
	 T[0]=0;
	 T[1]=1;
	 differ(S,T);
	 if ((X[1]!=T[0])||(X[2]!=T[1]))
	    error[1]+=1;
	 n=n+7;
	 count=count+1;
	 if ((numfac==2)&&(split==0)&&(psflag==0)) {
	    if ((yflag==0)||(yflag==1)) {
	       if (rflag!=sflag)
		  error[2]=8;
	       if ((rflag&0xf000)==0x2000)
		  error[2]=8;
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       quotient(T,T,ps);
	       bigprod(T[0],T[1],ps,X);
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       if ((X[1]!=T[0])||(X[2]!=T[1]))
		  error[2]=8;
	       }
	    else {
	       if ((rflag!=0xf333)&&(sflag!=0xf333)) {
		  if ((rflag&0xf000)!=0x2000)
		     error[2]=9;
		  if ((sflag&0xf000)!=0x2000)
		     error[2]=9;
		  }
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       quotient(T,T,ps);
	       bigprod(T[0],T[1],ps,X);
	       T[0]=0;
	       T[1]=1;
	       differ(W,T);
	       if ((X[1]!=T[0])||(X[2]!=T[1]))
		  error[2]=9;
	       }
	    }
	 }
askip: if (zflag==2)
          zflag=-1;
       zflag+=1;
       if (zflag==2)
          goto zloop;
}
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d count=%d error2=%d \n",error[0],error[1],error[2]);
fprintf(Outfp," count=%d \n",(n+1)/7);
for (i=0; i<(n+1)/7; i++)
   fprintf(Outfp," %#9x %#9x %#10x %#10x %#10x %#10x %#6x %#6x \n",
	   output[7*i],output[7*i+1],output[7*i+2],output[7*i+3],output[7*i+4],
	   output[7*i+5],(unsigned short)(output[7*i+6]>>16),
	   (unsigned short)(output[7*i+6]&0xffff));
if (error[2]!=0)
   printf(" error=%d \n",error[2]);
zskip:
fclose(Outfp);
return(0);
}