/*****************************************************************************/
/* */
/* 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 three distinct prime factors. */
/* */
/* The output is "a, b, (code0<<16)|code1, code2<<16, (("and" of code0, */
/* code1, and code2)<<16)|(((flag0<<2)|(flag1<<1)|flag2)<<4)|code". "code" */
/* is set to 0, 1, 2, or 3 if p divides a, b, a-b, or a+b respectively. */
/* "flag0" is set to 1 if the distinct prime factor of [(a**p+b**p)/(a+b)] */
/* is of form p**2*k+1, or 0 otherwise. "flag1" and "flag2" are similarly */
/* set. The corresponding bit in "code0" is set when p/2, 2p, 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), p**2*(a+b) */
/* is a pth power modulo a distinct prime factor of [(a**p+b**p)/(a+b)]. */
/* Simiarly, corresponding bits are set in "code1" and "code2". The three */
/* distinct prime factors of [(a**p+b**p)/(a+b)] are also output. */
/* */
/* 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.) */
/* */
/* Note: Not enough data has been collected for [(a**p+b**p)/(a+b)] to */
/* have three distinct prime factors of the form p**2*k+1. */
/* */
/*****************************************************************************/
#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 form p**2*k+1
// otherwise, factors can be of either form
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 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 error[];
unsigned int c,ps;
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=3500*7;
unsigned int d,e,temp,sumdif;
unsigned int i,j,k,l,m,oldk,oldoldk;
unsigned int flag,tflag,uflag,iters;
unsigned int aflag,bflag,cflag,dflag,limit;
int pflag,oflag;
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("out32c.dat","w");
ps=p*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;
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!=3)
goto askip;
if ((d/p)*p==d) {
dflag=0;
if (split==0) {
if ((d/2)*2!=d)
goto askip;
}
if (split==1) {
if ((d/2)*2==d)
goto askip;
}
}
if ((e/p)*p==e) {
dflag=1;
if (split==0) {
if ((e/2)*2!=e)
goto askip;
}
if (split==1) {
if ((e/2)*2==e)
goto askip;
}
}
if (((d+e)/p)*p==(d+e)) {
if (sumdif==1)
dflag=3;
else
dflag=2;
if (split==0) {
if (((d+e)/2)*2!=(d+e))
goto askip;
}
if (split==1) {
if (((d+e)/2)*2==(d+e))
goto askip;
}
}
if (((d-e)/p)*p==(d-e)) {
if (sumdif==1)
dflag=2;
else
dflag=3;
if (split==0) {
if (((d-e)/2)*2!=(d-e))
goto askip;
}
if (split==1) {
if (((d-e)/2)*2==(d-e))
goto askip;
}
}
/************************************/
/* 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 */
/************************************************/
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;
}
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))
flag=flag+4096;
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;
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;
}
if (((k-1)/ps)*ps==(k-1))
aflag=1;
else
aflag=0;
iters=iters+1;
if (iters==3)
break;
uflag=tflag;
tflag=flag;
cflag=bflag;
bflag=aflag;
oldoldk=oldk;
oldk=k;
}
if ((S[0]!=0)||(S[1]!=1)) {
error[0]=2;
goto bskip;
}
aflag=(aflag<<2)|(bflag<<1)|cflag;
oflag=uflag&tflag&flag;
flag=(flag<<16)|tflag;
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
if (n+6>outsiz) {
error[0]=6;
goto cskip;
}
output[n]=d;
output[n+1]=e;
output[n+2]=flag;
output[n+3]=(uflag<<16);
output[n+4]=(oflag<<16)|(aflag<<4)|dflag;
output[n+5]=(k<<16)|oldk;
output[n+6]=oldoldk;
n=n+7;
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," error=%d \n",error[0]);
fprintf(Outfp," count=%d \n",(n+1)/7);
for (i=0; i<(n+1)/7; i++)
fprintf(Outfp," %#10x %#10x %#10x %#10x %#10x %#10x %#6x \n",output[7*i],output[7*i+1],
output[7*i+2],output[7*i+3],output[7*i+4],output[7*i+5],output[7*i+6]);
fclose(Outfp);
return(0);
}