/*****************************************************************************/
/* */
/* FACTOR (a**p + b**p)/(a + b) */
/* 11/03/06 (dkc) */
/* */
/* This C program finds prime factors f of (a**p + b**p)/(a + b) where a */
/* and b are relatively prime integers and determines if q**((f-1)/p)=1 */
/* (mod f). All prime factors are output. A factor f must be duplicated. */
/* If p>3, whether p divides a when q divides a and q is a not a pth power */
/* modulo p**2 is determined (similarly for b). Whether p divides a-b or */
/* a+b when q divides a-b or a+b and q is not a pth power modulo p**2 is */
/* determined. */
/* */
/* The output is "a, b, [(a**p+b**p)/(a+b)] (two words), (n<<16)|code, */
/* factor1, factor1,...,factorn (possibly two words)" where "n" is the */
/* number of prime factors and "code" is set to p if p divides a+b, or */
/* 0 otherwise. If p does not divide a when q divides a, then an error */
/* is indicated ("error[1]" is set to a non-zero value). b is treated */
/* similarly. If p does not divide a-b or a+b when q divides a-b or a+b, */
/* then an error is indicated. */
/* */
/* Note: The user must ascertain that "base" is not a pth power modulo p^2. */
/* */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table12.h"
void product(unsigned int *a, unsigned int m, unsigned int count);
unsigned int lmbd(unsigned int mode, unsigned int a);
void sum(unsigned int *a, unsigned int *b);
void differ(unsigned int *a, unsigned int *b);
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 bigprod(unsigned int a, unsigned int b, unsigned int c, unsigned int *p);
int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 1000000.
// The maximum "dbeg" value for p=5 is about 5000.
// The maximum "dbeg" value for p=7 is about 500.
// The maximum "dbeg" value for p=11 is about 50.
//
unsigned int p=3; // input prime
unsigned int dbeg=1000000; // starting "a" value
unsigned int dend=1; // ending "a" value
//unsigned int stop=183;
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 base=12; // q value
// 1, 8 are pth powers modulo p**2 for p=3
// 1, 7, 18, 24 are pth powers modulo p**2 for p=5
// 1, 18, 19, 30, 31, 48 are pth powers modulo p**2 for p=7
// 1, 3, 9, 27, 40, 81, 94, 112, 118, 120 are pth powers
// modulo p**2 for p=11
extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int compos[];
extern unsigned int tmpsav;
extern unsigned int count;
extern unsigned int tcount;
extern unsigned int ccount;
extern unsigned int rflag;
unsigned int maxsiz=15600;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int outsiz=1999;
unsigned int cossiz=99;
unsigned int save[10]; // solutions array
unsigned int savsiz=9; // size of solutions array minus one
unsigned int d,e,a,b,temp;
unsigned int i,j,k,l,m;
unsigned int flag,dflag,bsave,limit;
unsigned int S[2],T[2],U[2],V[2],W[2],X[3];
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out23b.dat","w");
/*********************************/
/* extend prime look-up table */
/*********************************/
error[0]=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=2007; d<4000000; 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;
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
error[6]=0;
count=0;
tcount=0;
ccount=0;
rflag=0;
for (d=dbeg; d>=dend; d--) {
for (e=d-1; e>0; e--) {
// if (e!=stop) continue;
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;
if((d==(d/11)*11)&&(e==(e/11)*11)) continue;
if((d==(d/13)*13)&&(e==(e/13)*13)) continue;
if((d==(d/17)*17)&&(e==(e/17)*17)) continue;
if((d==(d/19)*19)&&(e==(e/19)*19)) 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 q divides d, e, d+e or d-e */
/******************************************/
bsave=0;
if (p!=3) {
if ((d/base)*base==d) {
bsave=d;
goto zskip;
}
if ((e/base)*base==e) {
bsave=e;
goto zskip;
}
}
if (((d+e)/base)*base==(d+e))
goto zskip;
if (((d-e)/base)*base!=(d-e))
continue;
/************************************/
/* compute (d**p + e**p)/(d + e) */
/************************************/
zskip:tcount=tcount+1;
dflag=0;
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];
}
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];
}
T[0]=X[1];
T[1]=X[2];
if (sumdif==1) {
sum(S, T);
temp=d+e;
if ((temp/p)*p==temp)
temp=temp*p;
quotient(T, S, temp);
}
else {
differ(S, T);
temp=d-e;
if ((temp/p)*p==temp)
temp=temp*p;
quotient(T, S, temp);
}
W[0]=S[0];
W[1]=S[1];
/************************************************/
/* look for prime factors using look-up table */
/************************************************/
if (S[0]==0) {
l = (33 - lmbd(1, S[1]))/2;
l = 1 << l;
}
else {
l = (65 - lmbd(1, S[0]))/2;
l = 1 << l;
}
k=0;
if (l>tmptab[tmpsiz-1]) {
flag=1;
k=tmpsiz-1;
}
else {
flag=0;
for (i=0; i<tmpsiz; i++) {
if (tmptab[i] < l) k=i;
else break;
}
}
l=k;
m=0;
for (i=0; i<=l; i++) {
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 (base!=1) {
bigresx(0, (k-1)/p, 0, k, U, base);
if ((U[0]!=0)||(U[1]!=1))
goto askip;
}
aloop: S[0]=V[0];
S[1]=V[1];
save[m]=k;
if (m < savsiz) m=m+1;
else {
error[0]=3;
goto bskip;
}
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])) {
dflag=1;
goto aloop;
}
}
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
if (dflag==0)
goto askip;
if ((S[0]!=0) || (S[1]!=1)) {
if (flag==1) {
if (S[0]==0) {
j = (33 - lmbd(1, S[1]))/2;
j = 1 << j;
}
else {
j = (65 - lmbd(1, S[0]))/2;
j = 1 << j;
}
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])) {
if (base!=1) {
bigresx(0, (i-1)/p, 0, i, U, base);
if ((U[0]!=0)||(U[1]!=1))
goto askip;
}
if (T[0]<=limit) { // largest prime in table is 1195153
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 {
if (ccount+3>cossiz) {
error[0]=4;
goto bskip;
}
compos[ccount]= (d<<16) | e;
compos[ccount+1] = S[0];
compos[ccount+2] = S[1];
compos[ccount+3] = i;
ccount = ccount+4;
goto askip;
}
}
}
}
cskip: if (base!=1) {
T[0]=0;
T[1]=1;
differ(S, T);
quotient(T, T, p);
bigresx(T[0], T[1], S[0], S[1], T, base);
}
else {
T[0]=0;
T[1]=1;
}
if ((T[0]==0)&&(T[1]==1)) {
if (n+m+6>outsiz) {
error[0]=6;
goto bskip;
}
if (m>0)
printf("d=%d, e=%d, m=%d \n",d,e,m+1);
output[n]=d;
output[n+1]=e;
output[n+2]=W[0];
output[n+3]=W[1];
if (sumdif==1) {
if (((d+e)/p)*p==(d+e)) k=p;
else k=0;
}
else {
if (((d-e)/p)*p==(d-e)) k=p;
else k=0;
}
output[n+4]=((int)(m+1) << 16) | (int)(k);
T[0]=S[0];
T[1]=S[1];
for (i=0; i<m; i++) {
bigprod(T[0], T[1], save[i], X);
T[0] = X[1];
T[1] = X[2];
output[n+i+5]=save[i];
}
output[n+m+5]=S[0];
output[n+m+6]=S[1];
if ((T[0]!=W[0]) || (T[1]!=W[1])) {
error[0]=7;
goto bskip;
}
n=n+m+7;
count=count+1;
if (bsave!=0) {
if ((bsave/p)*p!=bsave) {
error[1]+=1;
error[2]=d;
error[3]=e;
}
}
else {
if ((((d+e)/p)*p!=(d+e))&&(((d-e)/p)*p!=(d-e))) {
error[4]+=1;
error[5]=d;
error[6]=e;
}
}
}
else
goto askip;
}
else {
if (n+m+3>outsiz) {
error[0]=6;
goto bskip;
}
if (m>1)
printf("d=%d, e=%d, m=%d \n",d,e,m);
output[n]=d;
output[n+1]=e;
output[n+2]=W[0];
output[n+3]=W[1];
if (sumdif==1) {
if (((d+e)/p)*p==(d+e)) k=p;
else k=0;
}
else {
if (((d-e)/p)*p==(d-e)) k=p;
else k=0;
}
output[n+4]=((int)(m) << 16) | (int)(k);
S[0]=0;
S[1]=1;
for (i=0; i<m; i++) {
bigprod(S[0], S[1], save[i], X);
S[0] = X[1];
S[1] = X[2];
output[n+i+5]=save[i];
}
if ((S[0]!=W[0]) || (S[1]!=W[1])) {
error[0]=7;
goto bskip;
}
n=n+m+5;
count=count+1;
if (bsave!=0) {
if ((bsave/p)*p!=bsave) {
error[1]+=1;
error[2]=d;
error[3]=e;
}
}
else {
if ((((d+e)/p)*p!=(d+e))&&(((d-e)/p)*p!=(d-e))) {
error[4]+=1;
error[5]=d;
error[6]=e;
}
}
}
askip:rflag=1;
}
}
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," error4=%d asave=%d bsave=%d \n",error[4],error[5],
error[6]);
fprintf(Outfp," count=%d \n",n);
for (i=0; i<n; i++)
fprintf(Outfp," %#10x \n",output[i]);
fclose(Outfp);
if ((error[1]!=0)||(error[4]!=0))
printf(" error \n");
return(0);
}