/*****************************************************************************/
/* */
/* FACTOR (a**p+b**p)/(a+b) (when [(a**p+b**p)/(a+b)] is a pth power) */
/* 11/03/06 (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 */
/* prime factors of [(a**p+b**p)/(a+b)]. Use "table3b" (same as "table3a") */
/* for a<=1000000. Use "table4b" (same as "table4a" except for the range */
/* of a) for 1000000<a<=2000000. Use "table6b" (same as "table6a") for */
/* 2000000<a<=2500000. Use "table7b" (same as "table7a") for 2500000<a<= */
/* 3000000. Modify "insize" accordingly. */
/* */
/* Note: [(a**p+b**p)/(a+b)] can have only two distinct prime factors. */
/* */
/* The output is "a, b, (code0<<16)|code1, ("and" of code0 and code1)<<16| */
/* (count1<<2)|count0". The corresponding bit in "code0" is set when 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 distinct prime factor */
/* of [(a**p+b**p)/(a+b)]. Similarly, corresponding bits are set in */
/* "code1". "count1" is the number of distinct prime factors of */
/* [(a**p+b**p)/(a+b)] of the form p**2*k+1, and "count0" is the number */
/* of distinct prime factors of [(a**p+b**p)/(a+b)] not of the form */
/* p**2*k+1. The two distinct prime factors of [(a**p+b**p)/(a+b)] are */
/* also output. */
/* */
/* Additional output is "a, b, (flag0<<1)|flag1" when p is a pth power */
/* w.r.t. [(a**p+b**p)/(a+b)] and p**3 does not divide a, b, or a-b and */
/* p**4 does not divide a+b (the count should be zero when "split" is set */
/* to 0). "flag0" is set to 1 if the first distinct prime factor of */
/* [(a**p+b**p)/(a+b)] is of the form p**2*k+1, or 0 otherwise. "flag1" */
/* is similarly set for the second distinct prime factor of [(a**p+b**p)/ */
/* (a+b)]. */
/* */
/* Proposition (32) can be verified by examining the "code" when "psflag" */
/* is set to 1. (The "code" should also be examined when "mixed" types of */
/* factors of [(a**p+b**p)/(a+b)] are allowed.) */
/* */
/* Proposition (33) can be verified by examining the "code" when "psflag" */
/* is set to 0 and "split" is set to 0. (The "code" should also be examined*/
/* when "mixed" types of factors of [(a**p+b**p)/(a+b)] are allowed.) */
/* */
/* Proposition (35) can be verified by examining the "code" when "psflag" */
/* is set to 0 and "split" is set to 0. (The "code" should also be examined*/
/* when "mixed" types of factors of [(a**p+b**p)/(a+b)] are allowed.) */
/* */
/* Proposition (34) can be verified by examining the "code" when "psflag" */
/* is set to 0 and "split" is set to 1. (The "code" should also be examined*/
/* when "mixed" types of factors of [(a**p+b**p)/(a+b)] are allowed.) */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table3b.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);
void bigresx(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
unsigned int *e, unsigned int f);
int main ()
{
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 form p**2*k+1
// otherwise, factors can be of both types
unsigned int split=0; // 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 out=0; // if set, only output when p is a pth power w.r.t.
// (a^p+b^p)/(a+b)
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 tmpsav;
extern unsigned int rsave[];
extern unsigned int rcount;
extern unsigned int error[];
unsigned int c,ps,pc,pf;
unsigned int maxsiz=15000;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int insiz=1038; // table7b
//unsigned int insiz=1068; // table6b
//unsigned int insiz=2534; // table4b
//unsigned int insiz=4284; // table3b
unsigned int outsiz=5000*5;
unsigned int d,e,temp,limit;
unsigned int i,j,k,l,m,oldk;
unsigned int flag,rflag,sflag,tflag,uflag,vflag,iters,sumdif,aflag;
int pflag,qflag,count0,count1;
unsigned int S[2],T[2],U[2],V[2],X[3],Y[4],Z[4],Up[2];
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out32b.dat","w");
ps=p*p;
pc=ps*p;
pf=pc*p;
/*********************************/
/* 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]);
limit=(tmptab[tmpsiz-1])>>16;
limit=limit*limit;
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
pflag=0;
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
rcount=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!=2)
goto askip;
//
// check for "split" 2 and p
//
uflag=2;
if ((d/2)*2==d) {
if (sumdif!=0) {
if (((d+e)/p)*p==(d+e))
uflag=0;
}
else {
if (((d-e)/p)*p==(d-e))
uflag=0;
}
}
if ((e/2)*2==e) {
if (sumdif!=0) {
if (((d+e)/p)*p==(d+e))
uflag=1;
}
else {
if (((d-e)/p)*p==(d-e))
uflag=1;
}
}
if ((split==0)&&(uflag!=2))
goto askip;
if ((split==1)&&(uflag==2))
goto askip;
//
// check if p divides a, b, a-b, or a+b
//
if ((d/p)*p==d)
tflag=0;
if ((e/p)*p==e)
tflag=1;
if (((d+e)/p)*p==(d+e)) {
if (sumdif!=0)
tflag=3;
else
tflag=2;
}
if (((d-e)/p)*p==(d-e)) {
if (sumdif!=0)
tflag=2;
else
tflag=3;
}
//
// 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 */
/************************************************/
count0=0;
count1=0;
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 (psflag==1) {
if (((k-1)/ps)*ps!=(k-1))
goto askip;
}
if (psflag==0) {
if (((k-1)/ps)*ps==(k-1))
goto askip;
}
rflag=0;
if (((k-1)/ps)*ps==(k-1)) {
rflag=1;
count1+=1;
}
else
count0+=1;
flag=0;
bigresx(0, (k-1)/p, 0, k, U, d);
if ((U[0]==0)&&(U[1]==1))
flag=flag+2048;
bigresx(0, (k-1)/p, 0, k, U, e);
if ((U[0]==0)&&(U[1]==1))
flag=flag+1024;
bigresx(0, (k-1)/p, 0, k, U, d-e);
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+512;
else
flag=flag+256;
}
bigresx(0, (k-1)/p, 0, k, U, d+e);
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+256;
else
flag=flag+512;
}
bigresx(0, (k-1)/p, 0, k, U, p*d);
if ((U[0]==0)&&(U[1]==1))
flag=flag+128;
bigresx(0, (k-1)/p, 0, k, U, p*e);
if ((U[0]==0)&&(U[1]==1))
flag=flag+64;
bigresx(0, (k-1)/p, 0, k, U, p*(d-e));
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+32;
else
flag=flag+16;
}
bigresx(0, (k-1)/p, 0, k, U, p*(d+e));
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+16;
else
flag=flag+32;
}
bigresx(0, (k-1)/p, 0, k, U, p*p*d);
if ((U[0]==0)&&(U[1]==1))
flag=flag+8;
bigresx(0, (k-1)/p, 0, k, U, p*p*e);
if ((U[0]==0)&&(U[1]==1))
flag=flag+4;
bigresx(0, (k-1)/p, 0, k, U, p*p*(d-e));
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+2;
else
flag=flag+1;
}
bigresx(0, (k-1)/p, 0, k, U, p*p*(d+e));
if ((U[0]==0)&&(U[1]==1)) {
if (pflag!=2)
flag=flag+1;
else
flag=flag+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 (rflag!=0) {
if (flag!=0xfff)
error[1]=10;
}
else {
if (uflag==2) {
if (tflag==0) {
if (flag!=0x888)
error[1]=11;
}
if (tflag==1) {
if (flag!=0x444)
error[1]=11;
}
if (tflag>1) {
if (flag!=0x333)
error[1]=11;
}
}
}
flag=flag+4096;
}
else {
if (rflag!=0) {
if (tflag<3) {
if (flag!=0x10e)
error[1]=12;
}
else {
if (flag!=0xe10)
error[1]=12;
}
}
else {
if (uflag==2) {
if (tflag==0) {
if ((flag!=0x429)&&(flag!=0x258))
error[1]=13;
}
if (tflag==1) {
if ((flag!=0x825)&&(flag!=0x294))
error[1]=13;
}
if (tflag==2) {
if ((flag!=0x942)&&(flag!=0x582))
error[1]=13;
}
if (tflag==3) {
if ((flag!=0x294)&&(flag!=0x258))
error[1]=13;
}
}
else {
if (uflag==0) {
if ((flag!=0x942)&&(flag!=0x825))
error[1]=14;
}
if (uflag==1) {
if ((flag!=0x582)&&(flag!=0x429))
error[1]=14;
}
}
}
}
bigresx(0, (k-1)/p, 0, k, U, 2);
if ((Up[0]==U[0])&&(Up[1]==U[1]))
flag=flag+32768;
if ((U[0]==0)&&(U[1]==1))
flag=flag+8192;
bigresx(0, (k-1)/p, 0, k, U, 2*p);
if ((U[0]==0)&&(U[1]==1))
flag=flag+16384;
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]=1;
goto bskip;
}
iters=iters+1;
if (iters==2)
break;
sflag=flag;
vflag=rflag;
oldk=k;
}
if ((S[0]!=0)||(S[1]!=1)) {
error[0]=2;
goto bskip;
}
aflag=flag&sflag;
flag=(flag<<16)|sflag;
vflag=(vflag<<1)|rflag;
if ((aflag&0x1000)!=0) {
if (vflag==3) {
error[1]=15;
printf("warning: d=%d, e=%d, f=%d %d \n",d,e,k,oldk);
}
if (qflag==0) {
if (rcount<500) {
rsave[3*rcount]=d;
rsave[3*rcount+1]=e;
rsave[3*rcount+2]=vflag;
}
rcount+=1;
}
}
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
if (n+4>outsiz) {
error[0]=6;
goto cskip;
}
if ((out==0)||((aflag&0x1000)!=0)) {
output[n]=d;
output[n+1]=e;
output[n+2]=flag;
output[n+3]=(aflag<<16)|(count1<<2)|count0;
output[n+4]=(k<<16)|oldk;
n=n+5;
}
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)/5);
for (i=0; i<(n+1)/5; i++)
fprintf(Outfp," %#10x %#10x %#10x %#10x %#10x \n",output[5*i],output[5*i+1],
output[5*i+2],output[5*i+3],output[5*i+4]);
fprintf(Outfp,"\n");
fprintf(Outfp," p is a pth power w.r.t. [(a**p+b**p)/(a+b)] \n");
fprintf(Outfp," p**3 does not divide a, b, or a-b and p**4 does not divide a+b \n");
if (rcount<500)
n=rcount;
else
n=500;
fprintf(Outfp," count=%d \n",n);
for (i=0; i<n; i++) {
fprintf(Outfp," %#10x %#10x %#6x \n",rsave[3*i],rsave[3*i+1],rsave[3*i+2]);
}
fclose(Outfp);
if (error[1]!=0)
printf(" warning=%d \n",error[1]);
return(0);
}