/*****************************************************************************/
/* */
/* 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 one distinct prime factor. */
/* */
/* The output is "a, b, code, k". The corresponding bit in "code" 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 prime factor */
/* of [(a**p+b**p)/(a+b)]. "k" is the distinct prime factor of [(a^p+b^p)/ */
/* (a+b)]. */
/* */
/* Additional output is "a, b, m" where "m" is the number of prime factors */
/* of [(a**p+b**p)/(a+b)] when p is a pth power w.r.t. [(a**p+b**p)/(a+b)] */
/* and m is greater than 3. More output is "a, b, m" 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). More output is "a, b, m" when p is not a pth power w.r.t. */
/* and p**3 divides a, b, or a-b or p**4 divides a+b ("m" should be 9). */
/* The rest of Proposition (29) can be verified by examining the output */
/* "a, b, m" 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; p**3 should */
/* divide a', b', or a'-b' or p**4 should divide a'+b' for some */
/* representation [((a')**p+(b')**p)/(a'+b')] of T**p. */
/* */
/* 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 "table7b.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 not of the form p**2*k+1
// otherwise, factors can be of both forms
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 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 qsave[];
extern unsigned int rsave[];
extern unsigned int ssave[];
extern unsigned int tmpsav;
extern unsigned int count;
extern unsigned int qcount;
extern unsigned int rcount;
extern unsigned int scount;
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*3;
unsigned int d,e,temp;
unsigned int i,j,k,l,m;
unsigned int flag,sumdif,limit;
int pflag,qflag,rflag,tflag,uflag,t;
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("out32a.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) */
/***********************************/
count=0;
pflag=0;
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
qcount=0;
rcount=0;
scount=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!=1)
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;
t=(int)(d-2*e);
if (t<0)
t=-t;
if (((unsigned int)t/pc)*pc!=(unsigned int)t)
error[1]=88;
if (((2*d-e)/ps)*ps==(2*d-e))
error[1]=89;
}
}
else {
if (((d-e)/p)*p==(d-e)) {
uflag=0;
if (((d+2*e)/pc)*pc!=(d+2*e))
error[1]=88;
if (((2*d+e)/ps)*ps==(2*d+e))
error[1]=89;
}
}
}
if ((e/2)*2==e) {
if (sumdif!=0) {
if (((d+e)/p)*p==(d+e)) {
uflag=1;
if (((2*d-e)/pc)*pc!=(2*d-e))
error[1]=88;
if (((d-2*e)/ps)*ps==(d-2*e))
error[1]=89;
}
}
else {
if (((d-e)/p)*p==(d-e)) {
uflag=1;
if (((2*d+e)/pc)*pc!=(2*d+e))
error[1]=88;
if (((d+2*e)/ps)*ps==(d+2*e))
error[1]=89;
}
}
}
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;
if ((uflag==2)&&((2*d-e)/ps)*ps==(2*d-e))
error[1]=87;
t=(int)(d-2*e);
if (t<0)
t=-t;
if ((uflag==2)&&((unsigned int)t/ps)*ps==(unsigned int)t)
error[1]=87;
}
else
tflag=2;
}
if (((d-e)/p)*p==(d-e)) {
if (sumdif!=0)
tflag=2;
else {
tflag=3;
if ((uflag==2)&&((2*d+e)/ps)*ps==(2*d+e))
error[1]=87;
if ((uflag==2)&&((d+2*e)/ps)*ps==(d+2*e))
error[1]=87;
}
}
//
// 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 */
/************************************************/
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;
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]=9;
goto bskip;
}
else
break;
}
if ((m/3)*3!=m) {
error[0]=9;
goto bskip;
}
if ((S[0]!=0)||(S[1]!=1)) {
error[0]=9;
goto bskip;
}
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;
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;
if (m>3) {
if (rcount<500) {
rsave[3*rcount]=d;
rsave[3*rcount+1]=e;
rsave[3*rcount+2]=m;
}
rcount+=1;
}
if (qflag==0) {
if (scount<500) {
ssave[4*scount]=d;
ssave[4*scount+1]=e;
ssave[4*scount+2]=k;
ssave[4*scount+3]=m;
}
scount+=1;
}
}
else {
if (rflag!=0) {
if (tflag<3) {
if (flag!=0x10e)
error[1]=10;
}
else {
if (flag!=0xe10)
error[1]=10;
}
}
else {
if (uflag==2) {
if (tflag==0) {
if ((flag!=0x429)&&(flag!=0x258))
error[1]=12;
}
if (tflag==1) {
if ((flag!=0x825)&&(flag!=0x294))
error[1]=12;
}
if (tflag==2) {
if ((flag!=0x942)&&(flag!=0x582))
error[1]=12;
}
if (tflag==3) {
if ((flag!=0x294)&&(flag!=0x258))
error[1]=12;
}
}
else {
if (uflag==0) {
if ((flag!=0x942)&&(flag!=0x825))
error[1]=13;
}
if (uflag==1) {
if ((flag!=0x582)&&(flag!=0x429))
error[1]=13;
}
}
}
if (qflag==1) {
if (qcount<500) {
qsave[3*qcount]=d;
qsave[3*qcount+1]=e;
qsave[3*qcount+2]=m;
}
qcount+=1;
}
}
bigresx(0, (k-1)/p, 0, k, U, 2);
if ((U[0]==Up[0])&&(U[1]==Up[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;
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
if (n+2>outsiz) {
error[0]=6;
goto bskip;
}
if ((out==0)||((flag&0x1000)!=0)) {
output[n]=d;
output[n+1]=e;
output[n+2]=(flag<<16)|k;
n=n+3;
count+=1;
}
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 asave=%d bsave=%d csave=%d \n",
error[0],error[1],error[2],error[3],error[4]);
fprintf(Outfp," count=%d \n",(n+1)/3);
for (i=0; i<(n+1)/3; i++)
fprintf(Outfp," %#10x, %#10x, %#10x, %#10x, \n",output[3*i],output[3*i+1],
output[3*i+2]>>16,output[3*i+2]&0xffff);
fprintf(Outfp,"\n");
fprintf(Outfp," p is a pth power w.r.t. [(a**p+b**p)/(a+b)] \n");
fprintf(Outfp," [(a**p+b**p)/(a+b)] has more than 3 prime factors \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 %#10x \n",rsave[3*i],rsave[3*i+1],
rsave[3*i+2]);
}
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 (scount<500)
n=scount;
else
n=500;
fprintf(Outfp," count=%d \n",n);
for (i=0; i<n; i++) {
fprintf(Outfp," %#10x, %#10x, %#10x, %#10x, \n",ssave[4*i],ssave[4*i+1],
ssave[4*i+2],ssave[4*i+3]);
}
fprintf(Outfp,"\n");
fprintf(Outfp," p is not a pth power w.r.t. [(a**p+b**p)/(a+b)] \n");
fprintf(Outfp," p**3 divides a, b, or a-b, or p**4 divides a+b \n");
if (qcount<500)
n=qcount;
else
n=500;
fprintf(Outfp," count=%d \n",n);
for (i=0; i<n; i++) {
fprintf(Outfp," %#10x %#10x %#10x \n",qsave[3*i],qsave[3*i+1],
qsave[3*i+2]);
}
fclose(Outfp);
if (error[1]!=0)
printf(" errors=%#010x %d \n",error[0],error[1]);
return(0);
}