/*****************************************************************************/
/*									     */
/*  FACTOR (a**p+b**p)/(a+b)						     */
/*  11/03/06 (dkc)							     */
/*									     */
/*  This C program finds a and b such that (a**p + b**p)/(a + b) is a	     */
/*  square or p times a square.  p is set to 3.  p must divide a, b, or a-b  */
/*  or p**2 must divide a+b.  If p divides a, then q must divide a.  If p    */
/*  divides b, then q must divide b.  If p divides a-b or a+b, then q must   */
/*  divide a-b or a+b.	Whether q is a pth power modulus the square root of  */
/*  [(a**p+b**p)/(a+b)] is determined.	p should not divide q.		     */
/*									     */
/*  The output is "a, b".  If q is not a pth power modulus the square root   */
/*  of [(a**p+b**p)/(a+b)], then an error is indicated ("error[1]" is set    */
/*  to a non-zero value).						     */
/*									     */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table0a.h"
#include "table92.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
void dummy(unsigned int a, unsigned int b, unsigned int c);
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);
int main ()
{
unsigned int p=3;	 // input prime
unsigned int base=5;	 // q value

extern unsigned short table3[];
extern unsigned int input[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int insize;
extern unsigned int count;
extern unsigned int sumdif;
unsigned int t3size=2556;
unsigned int outsiz=1999;
unsigned int n=0;
unsigned int d,e,temp,save,ps,pb;
unsigned int h,i,j,k,l,lp,m;
unsigned int S[2],T[2],U[2],V[2],X[3];
double sqrt2,halfsr2;
FILE *Outfp;
Outfp = fopen("out22a3.dat","w");
sqrt2=1.4142135;
halfsr2=sqrt2*((double)(0.5));
/***********************************/
/*  factor (d**p + e**p)/(d + e)   */
/***********************************/
ps=p*p;
pb=p*base;
error[0]=0;	// clear error array
error[1]=0;
error[2]=0;
error[3]=0;
count=0;
   for (h=0; h<insize; h++) {
      d=input[2*h];
      e=input[2*h+1];
/******************************************/
/*  check if p divides d, e, d+e or d-e   */
/******************************************/
      if ((d/pb)*pb==d)
	 goto zskip;
      if ((e/pb)*pb==e)
	 goto zskip;
      if (sumdif==1) {
	 if (((d+e)/p)*p==(d+e)) {
	    if (((d+e)/ps)*ps!=(d+e))
	       continue;
	    }
	 }
      else {
	 if (((d-e)/p)*p==(d-e)) {
	    if (((d-e)/ps)*ps!=(d-e))
	       continue;
	    }
	 }
      if (((d+e)/p)*p==(d+e)) {
	 if (((d+e)/base)*base==(d+e))
	    goto zskip;
	 if (((d-e)/base)*base==(d-e))
	    goto zskip;
	 else
	    continue;
	 }
      if (((d-e)/p)*p==(d-e)) {
	 if (((d-e)/base)*base==(d-e))
	    goto zskip;
	 if (((d+e)/base)*base==(d+e))
	    goto zskip;
	 else
	    continue;
	 }
      else
	 continue;
/************************************/
/*  compute (d**p + e**p)/(d + e)   */
/************************************/
zskip:S[0]=0;
      S[1]=d;
      for (i=0; i<p-1; i++) {
	 bigprod(S[0], S[1], d, X);
	 S[0]=X[1];
	 S[1]=X[2];
	 }
      T[0]=0;
      T[1]=e;
      for (i=0; i<p-1; i++) {
	 bigprod(T[0], T[1], e, X);
	 T[0]=X[1];
	 T[1]=X[2];
	 }
      if (sumdif==1) {
	 sum(S, T);
	 temp=d+e;
	 if (((d+e)/p)*p==(d+e))
	    temp=temp*p;
	 quotient(T, S, temp);
	 }
      else {
	 differ(S, T);
	 temp=d-e;
	 if (((d-e)/p)*p==(d-e))
	    temp=temp*p;
	 quotient(T, S, temp);
	 }
/************************************************/
/*  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==0)
	 lp=(int)(((double)(l))*halfsr2);
      if (j==1) {
	 lp=l;
	 l=(int)(((double)(l))*sqrt2);
	 }
      lp=lp-1;
      l=l+1;
      if (l>table3[t3size-1]) {
	 error[0]=5;
	 goto bskip;
	 }
      else {
	 j=0;
	 for (i=0; i<t3size; i++) {
	    if (table3[i] < lp) j=i;
	    else break;
	    }
	 k=j;
	 for (i=j; i<t3size; i++) {
	    if (table3[i] < l) k=i;
	    else break;
	    }
	 }
      for (i=j; i<=k; i++) {
	 m=0;
	 l = table3[i];
	 quotient(S, V, l);
	 bigprod(V[0], V[1], l, X);
	 if ((S[0]!=X[1]) || (S[1]!=X[2])) continue;
	 save=l;
aloop:	 S[0]=V[0];
	 S[1]=V[1];
	 m=m+1;
	 quotient(S, V, l);
	 bigprod(V[0], V[1], l, X);
	 if ((S[0]==X[1]) && (S[1]==X[2])) goto aloop;
	 if (m!=2)
	    goto askip;
	 else
	    break;
	 }
      if (m!=2) continue;
      if ((S[0]!=0) || (S[1]!=1)) continue;
      bigresx(0, (save-1)/p, 0, save, U, base);
      if ((U[0]!=0)||(U[1]!=1)) {
	 error[1]=1;
	 error[2]=d;
	 error[3]=e;
	 goto bskip;
	 }
      if (n+1>outsiz) {
	 error[0]=6;
	 goto bskip;
	 }
      output[n]=d;
      output[n+1]=e;
      n=n+2;
      count=count+1;
askip:dummy(d, e, 0);
      }
bskip:
output[n]=-1;
fprintf(Outfp," error0=%d error1=%d asave=%d bsave=%d \n",error[0],
	error[1],error[2],error[3]);
fprintf(Outfp," count=%d \n",(n+1)/2);
for (i=0; i<(n+1)/2; i++)
   fprintf(Outfp," %#10x %#10x \n",output[2*i],output[2*i+1]);
fclose(Outfp);
if (error[1]!=0)
   printf(" error \n");
return(0);
}