/*****************************************************************************/
/* */
/* 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 */
/* fourth power or p times a fourth power. p is set to 3. p must divide */
/* a, b, a-b or a+b. If p divides a, q must divde a. If p divides b, q */
/* must divide b. If p divides a-b or a+b, q must divide a-b or a+b. */
/* Whether q is a pth power modulus the fourth 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 fourth 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 "table9.h"
#include "table94.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=7; // q value
extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int input[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int tmpsav;
extern unsigned int insize;
extern unsigned int count;
extern unsigned int sumdif;
unsigned int tsize=172;
unsigned int outsiz=1999;
unsigned int n=0;
unsigned int d,e,temp,save,flag,tmpsiz;
unsigned int h,i,j,k,l,lp,m,ps,pb;
unsigned int S[2],T[2],U[2],V[2],X[3];
double fourth2,fourth4,fourth8,half8;
FILE *Outfp;
Outfp = fopen("out22b3.dat","w");
fourth2=1.189207;
fourth4=1.4142135;
fourth8=1.6817928;
half8=0.8408964;
/*********************************/
/* 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;
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/4)*4;
l=l/4;
l = 1 << l;
if (j==0)
lp=(int)(((double)(l))*half8);
if (j==1) {
lp=l;
l=(int)(((double)(l))*fourth2);
}
if (j==2){
lp=(int)(((double)(l))*fourth2);
l=(int)(((double)(l))*fourth4);
}
if (j==3){
lp=(int)(((double)(l))*fourth4);
l=(int)(((double)(l))*fourth8);
}
lp=lp-1;
l=l+1;
if (l>tmptab[tmpsiz-1]) {
error[0]=5;
goto bskip;
}
else {
j=0;
for (i=0; i<tmpsiz; i++) {
if (tmptab[i] < lp) j=i;
else break;
}
k=j;
for (i=j; i<tmpsiz; i++) {
if (tmptab[i] < l) k=i;
else break;
}
}
for (i=j; i<=k; i++) {
m=0;
l = tmptab[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!=4)
goto askip;
else
break;
}
if (m!=4) 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);
}