/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/29/06 (dkc)							     */
/*									     */
/*  This C program determines if (a**p+b**p)/(a+b) is prime or p times a     */
/*  prime and of the form p**2*k+1.  Also, p**2 must divide a, b, a-b or     */
/*  a+b if p>3, or p**2 must divide a, b, a-b or p**3 must divide a+b when   */
/*  p=3.  Whether a, b, a-b, a+b, and p are pth powers modulus [(a**p+b**p)/ */
/*  (a+b)] is determined.						     */
/*									     */
/*  The output is "(a<<16)|b".  If a is not a pth power modulus [(a**p+b**p)/*/
/*  (a+b), "error[1]" is incremented.  Similarly, "error[2]", "error[3]",    */
/*  "error[4]", and "error[5]" are incremented if b, a-b, a+b, and p are     */
/*  not pth powers modulus [(a**p+b**p)/(a+b)]				     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table9.h"
void hugeprod(unsigned int a0, unsigned int a2, unsigned int a4,
	      unsigned int a6, unsigned int *product, unsigned int m0);
void dummy(unsigned int a, unsigned int b, unsigned int c);
unsigned int lmbd(unsigned int mode, unsigned int a);
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);
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 a, unsigned int b, unsigned int c, unsigned int d,
	     unsigned int *e, unsigned int f, unsigned int g);
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=150;	   // starting "a" value
unsigned int dend=1;	     // ending "a" value
//unsigned int stop=0x333;
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 base=2;  // base
unsigned int bflag=1; // if set, a, b, a-b, a+b, and p must be pth powers
		      // modulus [(a**p+b**p)/(a+b)], otherwise "base" must
		      // be a pth power modulus [(a**p+b**p)/(a+b)]

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 tsize=172;
unsigned int tmpsiz;
unsigned int outsiz=1000;
unsigned int d,e,a,b,temp;
unsigned int i,j,k,l,m;
unsigned int flag;
unsigned int S[2],T[2],U[2],X[3],Y[4],Z[4];
unsigned int ps,pc;
unsigned int n=0;
double sqrt2=1.4142135;
FILE *Outfp;
Outfp = fopen("out20.dat","w");
/*********************************/
/*  extend prime look-up table	 */
/*********************************/
for (i=0; i<tsize; i++) tmptab[i] = (int)(table[i]);
tmpsiz=tsize;
for (d=1033; d<170000; 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;
	 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;
	 }
   }
tmpsav=tmpsiz;
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p;
pc=ps*p;
error[0]=0;	// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
count=0;
for (d=dbeg; d>=dend; d--) {
   for (e=d-1; e>0; e--) {
//    if (e!=stop) continue;
/*********************************************/
/*  check if p**2 divides d, e, d+e or d-e   */
/*********************************************/
      if ((d/ps)*ps==d)
	 goto zskip;
      if ((e/ps)*ps==e)
	 goto zskip;
      if (p!=3) {
	 if (((d+e)/ps)*ps==(d+e))
	    goto zskip;
	 if (((d-e)/ps)*ps!=(d-e))
	    continue;
	 }
      else {
	 if (sumdif==1) {
	    if (((d+e)/pc)*pc==(d+e))
	       goto zskip;
	    if (((d-e)/ps)*ps!=(d-e))
	       continue;
	    }
	 else {
	    if (((d-e)/pc)*pc==(d-e))
	       goto zskip;
	    if (((d+e)/ps)*ps!=(d+e))
	       continue;
	    }
	 }
/******************************************/
/*  check for common factors of d and e   */
/******************************************/
zskip: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;
/************************************/
/*  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++)
	 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);
	 }
      S[0]=Y[2];
      S[1]=Y[3];
/************************************************/
/*  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/2)*2;
      l=l/2;
      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;
	    }
	 }
      for (i=0; i<=k; i++) {
	 m = tmptab[i];
	 quotient(S, T, m);
	 bigprod(T[0], T[1], m, X);
	 if ((S[0]==X[1]) && (S[1]==X[2]))
	    goto askip;
	 }
/***********************************************/
/*  output prime factors satisfying criterion  */
/***********************************************/
      if (flag==1) {
	 for (i=tmptab[tmpsiz-1]; i<=l; i+=2) {
	    quotient(S, T, i);
	    bigprod(T[0], T[1], i, X);
	    if ((X[1]==S[0]) && (X[2]==S[1]))
	       goto askip;
	    }
	 }
      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 (bflag==0) {
	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))
	   goto askip;
	}
      else {
	flag=0;
	T[0]=0;
	T[1]=1;
	differ(S, T);
	quotient(T, T, p);
	bigresx(T[0], T[1], S[0], S[1], U, d);
	if ((U[0]!=0)||(U[1]!=1)) {
	   error[1]+=1;
	   flag=1;
	   }
	bigresx(T[0], T[1], S[0], S[1], U, e);
	if ((U[0]!=0)||(U[1]!=1)) {
	   error[2]+=1;
	   flag=1;
	   }
	bigresx(T[0], T[1], S[0], S[1], U, d-e);
	if ((U[0]!=0)||(U[1]!=1)) {
	   error[3]+=1;
	   flag=1;
	   }
	bigresx(T[0], T[1], S[0], S[1], U, d+e);
	if ((U[0]!=0)||(U[1]!=1)) {
	   error[4]+=1;
	   flag=1;
	   }
	bigresx(T[0], T[1], S[0], S[1], U, p);
	if ((U[0]!=0)||(U[1]!=1)) {
	   error[5]+=1;
	   flag=1;
	   }
	if (flag!=0)
	   goto askip;
	}
      if (n>outsiz) {
	 error[0]=6;
	 n=n-1;
	 }
      output[n]=((int)(d) << 16) | (int)(e);
      n=n+1;
      count=count+1;
askip:dummy(d,e,8);
      }
   }
bskip:
output[n]=-1;
fprintf(Outfp," error0=%d error1=%d error2=%d error3=%d error4=%d error5=%d \n",
		error[0],error[1],error[2],error[3],error[4],error[5]);
fprintf(Outfp," count=%d \n",n);
if (n!=0) {
   for (i=0; i<n-1; i++)
      fprintf(Outfp," %#10x \n",output[i]);
   }
fclose(Outfp);
if ((error[1]!=0)||(error[2]!=0)||(error[3]!=0)||(error[4]!=0)||(error[5]!=0))
   printf(" error \n");
return(0);
}