/*****************************************************************************/
/* */
/* FACTOR (a**p+b**p)/(a+b) */
/* 11/29/06 (dkc) */
/* */
/* This C program finds a and b such that [(a**p+b**p)/(a+b)] and */
/* [(a**p-b**p)/(a-b)] are primes of the form p**2*k+1 and p**2 divides */
/* a, b, a-b, or a+b if p>3, or p**2 divides a or b, or p**3 divides a-b or */
/* a+b if p=3. Whether 2 is a pth power modulo [(a**p+b**p)/(a+b)] */
/* or [(a**p-b**p)/(a-b)] is determined. Whether [(a**p+b**p)/(a+b)] is */
/* a pth power modulo [(a**p-b**p)/(a-b)] and whether [(a**p-b**p)/(a-b)] */
/* is a pth power modulo [(a**p+b**p)/(a+b)] is determined. */
/* */
/* The output is "(a<<16)|b, (flag0<<16)|flag1". "flag0" consists of two */
/* bits; if the most significant bit is set, 2 is a pth power modulo */
/* [(a**p-b**p)/(a-b)] and if the least significant bit is set, 2 is a pth */
/* power modulo [(a**p+b**p)/(a+b)]. "flag1" consists of two bits; if */
/* the most significant bit is set, [(a**p+b**p)/(a+b)] is a pth power */
/* modulo [(a**p-b**p)/(a-b)] and if the least significant bit is set, */
/* [(a**p-b**p)/(a-b)] is a pth power modulo [[(a**p+b**p)/(a+b)]. If */
/* "flag0" is not equal to "flag1", an error is indicated ("error[1]" is */
/* is set to a non-zero value). */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "table3b.h"
void hugeprod(unsigned int a0, unsigned int a2, unsigned int a4,
unsigned int a6, unsigned int *product, unsigned int m0);
void dummy(unsigned int a, unsigned int b, unsigned int c);
unsigned int lmbd(unsigned int mode, unsigned int a);
void sum(unsigned int *addend, unsigned int *augend);
void differ(unsigned int *minuend, unsigned int *subtrahend);
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 hugeres(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
unsigned int *e, unsigned int f, unsigned int g);
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);
int main ()
{
//
// Note: The maximum "dbeg" value for p=3 is about 65000.
// The maximum "dbeg" value for p=5 is about 65000.
// The maximum "dbeg" value for p=7 is about 1000.
// The maximum "dbeg" value for p=11 is about 50.
//
unsigned int p=3; // prime value
unsigned int dbeg=1000; // starting "a" value
unsigned int dend=1; // ending "a" value
unsigned int rflag=1; // if set, p**2 must divide d, e, d+e, or d-e
//unsigned int stop=0xf4;
extern unsigned short table[]; // prime look-up table
extern unsigned int tmptab[]; // extended prime look-up table
extern unsigned int output[]; // output array
extern unsigned int error[]; // error array
unsigned int tmpsav; // size of extended prime look-up table
extern unsigned int count; // output count
unsigned int tsize=1228; // size of prime look-up table
unsigned int tmpsiz; // size of extended prime look-up table
unsigned int outsiz=1000*2; // size of output array
unsigned int d,e;
unsigned int a,b,temp;
unsigned int i,j,k,l;
unsigned int flag,sumdif,ps,pc;
unsigned int tflag,sflag;
unsigned int P[4],Q[2],R[2],S[2],T[2],U[2],X[3],Y[4],Z[4];
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out16.dat","w");
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
/*********************************/
/* extend prime look-up table */
/*********************************/
//
// copy prime look-up table
//
for (i=0; i<tsize; i++) tmptab[i] = (unsigned int)(table[i]);
tmpsiz=tsize;
//
// check for small prime factors
//
for (d=10001; d<160000; d++) {
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;
//
// find size of largest possible prime factor
// find maximum index into prime look-up table
//
l = (int)(100.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; // maximum index needed
//
// check for prime factors using look-up table
//
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;
}
}
tmpsav=tmpsiz;
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
ps=p*p; // p**2
pc=ps*p; // p**3
count=0;
for (d=dbeg; d>=dend; d--) { // loop through d values
for (e=d-1; e>0; e--) { // loop through e values
// if (e!=stop) continue;
/*********************************************/
/* check if p**2 divides d, e, d+e or d-e */
/*********************************************/
if (rflag!=0) {
if ((d/ps)*ps==d) // check if p**2 divides d
goto zskip;
if ((e/ps)*ps==e) // check if p**2 divides e
goto zskip;
if (p!=3) {
if (((d+e)/ps)*ps==(d+e)) // check if p**2 divides d+e
goto zskip;
if (((d-e)/ps)*ps!=(d-e)) // check if p**2 divides d-e
continue;
}
else {
if (((d+e)/pc)*pc==(d+e)) // check if p**3 divides d+e
goto zskip;
if (((d-e)/pc)*pc!=(d-e)) // check if p**e divides d-e
continue;
}
}
/****************************************/
/* check for common factors of d and e */
/****************************************/
zskip: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;
/*******************************/
/* compute (d**p+e**p)/(d+e) */
/*******************************/
tflag=0;
sflag=0;
Y[0]=0;
Y[1]=0;
Y[2]=0;
Y[3]=d;
for (i=0; i<p-1; i++)
hugeprod(Y[0], Y[1], Y[2], Y[3], Y, d); // multiple-word product
Z[0]=0;
Z[1]=0;
Z[2]=0;
Z[3]=e;
for (i=0; i<p-1; i++)
hugeprod(Z[0], Z[1], Z[2], Z[3], Z, e); // multiple-word product
P[0]=Z[0];
P[1]=Z[1];
P[2]=Z[2];
P[3]=Z[3];
bigbigs(Y, Z); // multiple-word sum
bigbigd(Y, P); // multiple-word difference
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); // multiple-word quotient
R[0]=Y[2];
R[1]=Y[3];
temp=d-e;
if (((d-e)/p)*p==(d-e))
temp=temp*p;
bigbigq(P[0], P[1], P[2], P[3], Y, 0, temp); // multiple-word quotient
Q[0]=Y[2];
Q[1]=Y[3];
/***********************/
/* begin inner loop */
/***********************/
for (sumdif=0; sumdif<2; sumdif++) {
if (sumdif==0) { // select between q and r
S[0]=Q[0];
S[1]=Q[1];
}
else {
S[0]=R[0];
S[1]=R[1];
}
/*********************************************************/
/* find size of largest possible prime factor */
/* find maximum index into extended prime look-up table */
/*********************************************************/
if (S[0]==0) {
l = (33 - lmbd(1, S[1]))/2; // left-most bit detection
l = 1 << l;
}
else {
l = (65 - lmbd(1, S[0]))/2; // left-most bit detection
l = 1 << l;
}
k=0;
if (l>tmptab[tmpsiz-1]) {
k=tmpsiz-1;
}
else {
for (i=0; i<tmpsiz; i++) {
if (tmptab[i] < l) k=i;
else break;
}
}
l=k; // maximum index needed
/************************************************/
/* look for prime factors using look-up table */
/************************************************/
for (i=0; i<=l; i++) {
k = tmptab[i]; // load prime
quotient(S, T, k); // multiple-word quotient
bigprod(T[0], T[1], k, X); // multiple-word product
if ((S[0]==X[1]) && (S[1]==X[2])) goto askip; // check if factor
}
/***********************************************************/
/* continue looking for factors (not necessarily prime) */
/***********************************************************/
//
// find size of largest possible factor
//
if (S[0]==0) {
j = (33 - lmbd(1, S[1]))/2; // left-most bit detection
j = 1 << j;
}
else {
j = (65 - lmbd(1, S[0]))/2; // left-most bit detection
j = 1 << j;
}
//
// loop through potential factors
//
for (i=tmptab[tmpsiz-1]; i<j; i+=2) {
quotient(S, T, i); // multiple-word quotient
bigprod(T[0], T[1], i, X); // multiple-word product
if ((X[1]==S[0]) && (X[2]==S[1])) // check if factor
goto askip;
}
/****************************************************************************/
/* check if [(a**p+b**p)/(a+b)] or [(a**p-b**p)/(a-b)] is of form p**2*k+1 */
/****************************************************************************/
T[0]=0;
T[1]=1;
differ(S, T); // multiple-word difference
quotient(T, T, ps); // multiple-word quotient
bigprod(T[0], T[1], ps, X); // multiple-word product
T[0]=0;
T[1]=1;
differ(S, T); // multiple-word difference
if ((X[1]!=T[0]) || (X[2]!=T[1])) // check if p**2 is factor
goto askip;
/***********************************/
/* check if proposition satisfied */
/***********************************/
quotient(T, T, p); // (q-1)/p
hugeres(T[0], T[1], S[0], S[1], U, 0, 2); // 2**((q-1)/p) modulus q
if ((U[0]!=0)||(U[1]!=1)) // check if 2**((q-1)/p)==1(mod q)
tflag=tflag|1; // set flag
if (sumdif==0)
hugeres(T[0], T[1], S[0], S[1], U, R[0], R[1]); // least residue
else
hugeres(T[0], T[1], S[0], S[1], U, Q[0], Q[1]); // least residue
if ((U[0]!=0)||(U[1]!=1)) // check if r**((q-1)/p)=1(mod q)
sflag=sflag|1; // set flag
if (sumdif==0) {
tflag=tflag<<1; // left-shift flag for packing
sflag=sflag<<1; // left-shift flag for packing
}
else { // output results
if (n+1>outsiz) {
error[0]=6;
goto bskip;
}
output[n]=(d<<16) | e; // output d and e
output[n+1]=(tflag<<16)|sflag; // output packed flags
n=n+2;
count=count+1; // increment count
if (tflag!=sflag)
error[1]=8;
}
}
askip: dummy(d,e,0);
}
}
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d error1=%d \n",error[0],error[1]);
fprintf(Outfp," count=%d \n",(n+1)/2);
for (i=0; i<(n+1)/2; i++)
fprintf(Outfp," %#10x %#10x \n",output[2*i],output[2*i+1]);
fclose(Outfp);
if (error[1]!=0)
printf(" error \n");
return(0);
}