/*****************************************************************************/
/* */
/* FACTOR (a**p + b**p)/(a + b) */
/* 03/15/14 (dkc) */
/* */
/* This C program determines if p is a pth power with respect to */
/* (a**p + b**p)/(a+b). p must divide a, b, a-b or a+b. Whether p**2 */
/* divides a, b, a-b, or a+b is determined. Whether the prime factors of */
/* [(a^p+b^p)/(a+b)] are of the form p^2k+1 is determined. */
/* */
/* The output is "(a<<16)|b". If p**2 does not divide a, b, or a-b, then */
/* an error is indicated ("error[1]" is set to a non-zero value). If p>3 */
/* and p**2 does not divide a+b, then an error is indicated. If p=3 and */
/* p**3 does not divide a+b when p**2 divides a+b, then an error is */
/* indicated. */
/* */
/* Note: Prime [(a^p+b^p)/(a+b)] are sometimes not checked for "very large" */
/* [(a^p+b^p)/(a+b)] values. */
/* */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table12.h"
unsigned int lmbd(unsigned int mode, unsigned int a);
void bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
void bigbigd(unsigned int *a, unsigned int *b);
void differ(unsigned int *a, unsigned int *b);
void dummy(unsigned int a, unsigned int b, unsigned int c);
void bigbigs(unsigned int *addend, unsigned int *augend);
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);
void quotient(unsigned int *a, unsigned int *b, unsigned int c);
void newprod(unsigned int a, unsigned int b, unsigned int c, unsigned *d,
unsigned int e);
int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 2^32.
// The maximum "dbeg" value for p=5 is about 2^20.
// The maximum "dbeg" value for p=7 is about 2^14.
// The maximum "dbeg" value for p=11 is about 2^9.
//
unsigned int p=7; // input prime
unsigned int dbeg=1000; // starting "a" value
unsigned int dend=1; // ending "a" value
//unsigned int stop=57;
unsigned int sumdif=1; // select [(a**p+b**p)/(a+b)] if "sumdif" is non-
// zero, or [(a**p-b**p)/(a-b)] otherwise
unsigned int correct=0;
// There is a small probability that (a**p+b**p)/(a+b) is not
// completely factored if "correct" is not set to 1.
unsigned int out=0; // output flag
extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int tmpsav;
extern unsigned int count;
unsigned int maxsiz=50000;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int outsiz=1999;
unsigned int save[16]; // solutions array
unsigned int savsiz=15; // size of solutions array minus one
unsigned int d,e,a,b,temp,lcount,scount,tcount;
unsigned int i,j,k,l,m;
unsigned int flag,limit,dflag,xflag,pflag,qflag;
unsigned int S[4],T[4],U[4],V[4],W[4],X[4],Y[4],Z[4];
unsigned int ps,pc,minfac;
unsigned int n=0;
double sqrt2=1.4142135;
FILE *Outfp;
Outfp = fopen("out42.dat","w");
/*********************************/
/* extend prime look-up table */
/*********************************/
error[0]=0;
error[1]=0;
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=9975; d<(9971*9971); 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;
}
}
tmpsav=tmpsiz;
limit=(tmptab[tmpsiz-1])>>16;
limit=limit*limit;
printf("count=%d, prime=%d \n",tmpsav,tmptab[tmpsiz-1]);
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
minfac=0x7fffffff;
lcount=0;
scount=0;
tcount=0;
ps=p*p;
pc=ps*p;
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
error[6]=0;
error[7]=0;
error[8]=0;
error[9]=0;
count=0;
for (d=dbeg; d>=dend; d--) {
for (e=d-1; e>0; e--) {
// if (e!=stop) continue;
/******************************************/
/* check for common factors of d and e */
/******************************************/
if((d==(d/2)*2)&&(e==(e/2)*2)) continue;
if((d==(d/3)*3)&&(e==(e/3)*3)) continue;
if((d==(d/5)*5)&&(e==(e/5)*5)) continue;
if((d==(d/7)*7)&&(e==(e/7)*7)) continue;
/***********************/
/* Euclidean G.C.D. */
/***********************/
a=d;
b=e;
if (b>a) {
temp=a;
a=b;
b=temp;
}
loop: temp = a - (a/b)*b;
a=b;
b=temp;
if (b!=0) goto loop;
if (a!=1) continue;
/********************************************/
/* check if p**2 divide a, b, a+b or a-b */
/********************************************/
xflag=0;
if ((d/p)*p==d) {
xflag=d;
goto jump;
}
if ((e/p)*p==e) {
xflag=e;
goto jump;
}
if (sumdif!=0) {
if (((d-e)/p)*p==(d-e)) {
xflag=d-e;
goto jump;
}
if (p==3) {
if (((d+e)/ps)*ps!=(d+e))
continue;
else
xflag=d+e;
}
else {
if (((d+e)/p)*p!=(d+e))
continue;
else
xflag=d+e;
}
}
else {
if (((d+e)/p)*p==(d+e)) {
xflag=d+e;
goto jump;
}
if (p==3) {
if (((d-e)/ps)*ps!=(d-e))
continue;
else
xflag=d-e;
}
else {
if (((d-e)/p)*p!=(d-e))
continue;
else
xflag=d-e;
}
}
jump: dflag=0;
if ((d/ps)*ps==d) {
dflag=d;
goto zskip;
}
if ((e/ps)*ps==e) {
dflag=e;
goto zskip;
}
if (sumdif!=0) {
if (((d-e)/ps)*ps==(d-e)) {
dflag=d-e;
goto zskip;
}
if (p==3) {
if (((d+e)/pc)*pc==(d+e))
dflag=d+e;
}
else {
if (((d+e)/ps)*ps==(d+e))
dflag=d+e;
}
}
else {
if (((d+e)/ps)*ps==(d+e)) {
dflag=d+e;
goto zskip;
}
if (p==3) {
if (((d-e)/pc)*pc==(d-e))
dflag=d-e;
}
else {
if (((d-e)/ps)*ps==(d-e))
dflag=d-e;
}
}
/************************************/
/* compute (d**p + e**p)/(d + e) */
/************************************/
zskip:Y[0]=0;
Y[1]=0;
Y[2]=0;
Y[3]=d;
for (i=0; i<p-1; i++)
newprod(Y[1],Y[2],Y[3],Y,d);
Z[0]=0;
Z[1]=0;
Z[2]=0;
Z[3]=e;
for (i=0; i<p-1; i++)
newprod(Z[1],Z[2],Z[3],Z,e);
if (sumdif==1) {
bigbigs(Y, Z);
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);
}
else {
bigbigd(Y, Z);
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);
}
W[0]=Y[0];
W[1]=Y[1];
W[2]=Y[2];
W[3]=Y[3];
/************************************************/
/* look for prime factors using look-up table */
/************************************************/
if (Y[0]!=0) {
printf("S value too big \n");
goto zskip;
}
if (Y[1]!=0)
l=96-lmbd(1,Y[1]);
else {
if (Y[2]!=0)
l=64-lmbd(1,Y[2]);
else
l=32-lmbd(1,Y[3]);
}
j=l-(l/2)*2;
l=l/2;
if (l>28) {
flag=1;
k=tmpsiz-1;
goto ajump;
}
l = 1 << l;
if (j==1)
l=(int)(((double)(l))*sqrt2);
l=l+1;
flag=0;
if (l>tmptab[tmpsiz-1]) {
flag=1;
k=tmpsiz-1;
}
else {
k=0;
for (i=0; i<tmpsiz; i++) {
if (tmptab[i] < l) k=i;
else break;
}
}
ajump:m=0;
for (i=0; i<savsiz; i++)
save[i]=0;
for (i=0; i<=k; i++) {
l = tmptab[i];
bigbigq(Y[0],Y[1],Y[2],Y[3],T,0,l);
V[0]=T[0];
V[1]=T[1];
V[2]=T[2];
V[3]=T[3];
newprod(T[1],T[2],T[3],T,l);
if ((Y[0]!=T[0])||(Y[1]!=T[1])||(Y[2]!=T[2])||(Y[3]!=T[3]))
continue;
aloop: Y[0]=V[0];
Y[1]=V[1];
Y[2]=V[2];
Y[3]=V[3];
save[m]=l;
if (m < savsiz) m=m+1;
else {
error[0]=3;
goto bskip;
}
bigbigq(Y[0], Y[1], Y[2], Y[3], T, 0, l);
V[0]=T[0];
V[1]=T[1];
V[2]=T[2];
V[3]=T[3];
newprod(T[1],T[2],T[3],T,l);
if ((Y[0]==T[0])&&(Y[1]==T[1])&&(Y[2]==T[2])&&(Y[3]==T[3]))
goto aloop;
}
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
if ((Y[0]!=0)||(Y[1]!=0)||(Y[2]>0x3fffffff))
goto askip;
S[0]=Y[2];
S[1]=Y[3];
if ((S[0]!=0) || (S[1]!=1)) {
if ((flag==1) && (correct==1)) {
if (S[0]==0)
j = (32 - lmbd(1, S[1]));
else
j = (64 - lmbd(1, S[0]));
k=j-(j/2)*2;
j=j/2;
j = 1 << j;
if (k==1)
j=(int)(((double)(j))*sqrt2);
for (i=tmptab[tmpsiz-1]; i<j; i+=2*p) {
quotient(S, T, i);
bigprod(T[0], T[1], i, X);
if ((X[1]==S[0]) && (X[2]==S[1])) {
bigresx(0, (i-1)/p, 0, i, U, p);
if ((U[0]!=0)||(U[1]!=1)) {
goto askip;
}
if (T[0]<=limit) { // largest prime in table is 0x126f5f
S[0]=T[0]; // for p=7
S[1]=T[1];
save[m]=i;
if (m < savsiz) m=m+1;
else {
error[0]=3;
goto bskip;
}
goto cskip;
}
else {
error[0]=4;
goto bskip;
}
}
}
}
cskip: T[0]=0;
T[1]=1;
differ(S, T);
quotient(T, T, p);
bigresx(T[0], T[1], S[0], S[1], U, p);
if ((U[0]==0)&&(U[1]==1)) {
//
// check if of the form p^2k+1
//
qflag=1;
T[0]=0;
T[1]=1;
differ(S, T);
quotient(T, U, ps);
bigprod(U[0],U[1],ps,X);
if ((X[1]==T[0])&&(X[2]==T[1]))
qflag=0;
//
T[0]=0;
T[1]=0;
T[2]=S[0];
T[3]=S[1];
pflag=1;
for (i=0; i<m; i++) {
l=save[i];
bigresx(0, (l-1)/p, 0, l, U, p);
if ((U[0]!=0)||(U[1]!=1))
goto askip;
newprod(T[1],T[2],T[3],T,l);
if (((l-1)/ps)*ps==(l-1)) {
pflag=0;
if (l<minfac)
minfac=l;
}
}
if ((T[0]!=W[0])||(T[1]!=W[1])||(T[2]!=W[2])||(T[3]!=W[3])) {
error[0]=7;
goto bskip;
}
if (out!=0) {
if (n>outsiz) {
error[0]=6;
goto bskip;
}
output[n]=((int)(d) << 16) | (int)(e);
n=n+1;
}
if (m>0) {
count=count+1;
tcount=tcount+m+1;
printf("d=%d, e=%d, m=%d, f=%d %d %d %d %d %d %#010x %#010x \n",d,e,m+1,save[0],save[1],save[2],save[3],save[4],
save[6], S[0],S[1]);
fprintf(Outfp,"d=%d, e=%d, m=%d, f=%d %d %d %d %d %d %#010x %#010x \n",d,e,m+1,save[0],save[1],save[2],save[3],save[4],
save[6], S[0],S[1]);
if (pflag==0) {
printf("small p^2*k+1: d=%d, e=%d \n",d,e);
fprintf(Outfp,"small p^2*k+1: d=%d, e=%d \n",d,e);
scount=scount+1;
}
if (qflag==0) {
printf("large p^2*k+1: d=%d, e=%d \n",d,e);
fprintf(Outfp,"large p^2*k+1: d=%d, e=%d \n",d,e);
lcount=lcount+1;
}
if ((pflag==0)&&(qflag==0)) {
printf("warning: both flags set \n");
fprintf(Outfp,"warning: both flags set \n");
}
}
if (dflag!=xflag) {
error[1]=8;
printf("error: d=%d, e=%d, flags=%d %d \n",d,e,dflag,xflag);
goto bskip;
}
}
else {
goto askip;
}
}
else {
T[0]=0;
T[1]=0;
T[2]=0;
T[3]=1;
pflag=1;
for (i=0; i<m; i++) {
l=save[i];
bigresx(0, (l-1)/p, 0, l, U, p);
if ((U[0]!=0)||(U[1]!=1))
goto askip;
if ((m>1)&&(i!=(m-1))) {
if (((l-1)/ps)*ps==(l-1)) {
pflag=0;
if (l<minfac)
minfac=l;
}
}
newprod(T[1],T[2],T[3],T,l);
}
qflag=1;
if (m>1) {
l=save[m-1];
if (((l-1)/ps)*ps==(l-1))
qflag=0;
}
if ((T[0]!=W[0])||(T[1]!=W[1])||(T[2]!=W[2])||(T[3]!=W[3])) {
error[0]=7;
goto bskip;
}
if (out!=0) {
if (n>outsiz) {
error[0]=6;
goto bskip;
}
output[n]=((int)(d) << 16) | (int)(e);
n=n+1;
}
if (m>1) {
count=count+1;
tcount=tcount+m;
printf("d=%d, e=%d, m=%d, f=%d %d %d %d %d %d \n",d,e,m,save[0],save[1],save[2],save[3],save[4],
save[6]);
fprintf(Outfp,"d=%d, e=%d, m=%d, f=%d %d %d %d %d %d \n",d,e,m,save[0],save[1],save[2],save[3],save[4],
save[6]);
if (pflag==0) {
printf("small p^2*k+1: d=%d, e=%d \n",d,e);
fprintf(Outfp,"small p^2*k+1: d=%d, e=%d \n",d,e);
scount=scount+1;
}
if (qflag==0) {
printf("large p^2*k+1: d=%d, e=%d \n",d,e);
fprintf(Outfp,"large p^2*k+1: d=%d, e=%d \n",d,e);
lcount=lcount+1;
}
if ((pflag==0)&&(qflag==0)) {
printf("warning: both flags set \n");
fprintf(Outfp,"warning: both flags set \n");
}
}
if (dflag!=xflag) {
error[1]=8;
printf("error: d=%d, e=%d, flags=%d %d \n",d,e,dflag,xflag);
goto bskip;
}
}
askip:dummy(d,e,6);
}
}
bskip:
if (out!=0)
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d \n",error[0],error[1]);
fprintf(Outfp," count=%d \n",count);
printf(" count=%d \n",count);
if ((n!=0)&&(out!=0)) {
for (i=0; i<n-1; i++)
fprintf(Outfp," %#10x \n",output[i]);
}
if (error[1]!=0)
printf(" error=%d \n",error[1]);
printf("small=%d, large=%d, total=%d, minimum=%d \n",scount,lcount,tcount,minfac);
fprintf(Outfp,"small=%d, large=%d, total=%d, minimum=%d \n",scount,lcount,tcount,minfac);
fclose(Outfp);
return(0);
}