/*****************************************************************************/
/* */
/* FACTOR (a**p-b**p)/(a-b) (when [(a**p+b**p)/(a+b)] is a pth power) */
/* 11/03/06 (dkc) */
/* */
/* This C program determines if q is a pth power with respect to */
/* (a**p-b**p)/(a-b). Use "table4a.h" for a<=2000000 and "table6a.h" for */
/* for 2000000<a<=2500000. Use "table7a" for 2500000<a<=3000000. Modify */
/* "insize" accordingly. */
/* */
/* The output is "a, b, (number of prime factors<<16)|code". The code is */
/* set to 0, 1, 2, or 3 when p divides a, b, a-b, or a+b respectively. */
/* Setting "select" to 0, 1, 2, 3, 4, 5, or 6 selects a, b, a-b, a+b, p, */
/* 2, or 2*p respectively as the base (q). When "select" ranges from */
/* 0 to 3, setting "pflag" to 1 multplies the base by p and setting "pflag" */
/* to 2 multiplies the base by p*p. 2*p must divide a, b, a-b, or a+b */
/* when "split" is set to 0, and 2*p must not divide a, b, a-b, or a+b when */
/* "split" is set to 1. The prime factors of [(a**p-b**p)/(a-b)] must not */
/* be of the form p**2*k+1 when "psflag" is set to 0, and the prime factors */
/* of [(a**p-b**p)/(a-b)] must be of the form p**2*k+1 when "psflag" is set */
/* to 1. */
/* */
/* When "split" equals 0 and "select" equals 2 or 3, "pflag" should be set */
/* to 1. When "split" equals 1 and "select" equals 0 or 1, "pflag" should */
/* be set to 1. */
/* */
/* When "split" equals 0, "code" equals "select", and "base" is not a pth */
/* power w.r.t [(a**p-b**p)/(a-b)], then an error is indicated ("error[0]" */
/* is set to 8). This confirms Proposition (12). */
/* */
/* When "split" equals 0, "select" equals 5, "code" equals 2 or 3, and */
/* "base" is not a pth power w.r.t. [(a**p-b**p)/(a-b)], then an error is */
/* is indicated ("error[0]" is set to 9). This confirms Proposition (13). */
/* */
/* When "split" equals 1 and "select" equals 5, the output count is 0. */
/* This confirms part of Proposition (3). */
/* */
/* When "split" equals 1 and "select" equals 4, the output count is 0. */
/* This confirms the rest of Proposition (3). */
/* */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table7a.h"
void differ(unsigned int *a, unsigned int *b);
void bigbigs(unsigned int *a, unsigned int *b);
void bigbigd(unsigned int *a, unsigned int *b);
void bigbigq(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int *quotient, unsigned int d2,
unsigned int d3);
unsigned int lmbd(unsigned int mode, unsigned int a);
void dummy(unsigned int a, unsigned int b, unsigned int c);
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 bigresx(unsigned int a, unsigned int b, unsigned int c, unsigned int d,
unsigned int *e, unsigned int f);
int main ()
{
unsigned int p=3; // input prime
unsigned int select=4; // select d, e, d-e, d+e, p, 2, or 2*p
unsigned int pflag=0; // if set to 1, base = base*p for "select" = 0 to 3
// if set to 2, base = base*p*p for "select"=0 to 3
unsigned int split=1; // if set to 0, 2*p must divide a, b, a-b, or a+b
// if set to 1, 2*p must not divide a, b, a-b, or a+b
// otherwise, no restrictions
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 of the form p**2*k+1
// otherwise, factors can be of mixed types
unsigned int qflag=4; // if set to 0, p must divide a+b
// if set to 1, p must divide a-b
// if set to 2, p must divide a
// if set to 3, p must divide b
// otherwise, no restrictions
unsigned int prmflg=0; // if set to 1, [(a**p-b**p)/(a-b)] must be prime
unsigned int offset=0; // offset into the input look-up table (must be
// even and less than "insize"). Used to reduce
// execution time.
extern unsigned int input[];
extern unsigned short table[];
extern unsigned int tmptab[];
extern unsigned int output[];
extern unsigned int terror[];
extern unsigned int count;
unsigned int insize=866; // table7a
//unsigned int insize=1068; // table6a
//unsigned int insize=6818; // table4a
unsigned int maxsiz=200000;
unsigned int tsize=1228;
unsigned int tmpsiz;
unsigned int outsiz=10000*3;
unsigned int save[20]; // solutions array
unsigned int savsiz=19; // size of solutions array minus one
unsigned int d,e,c;
unsigned int h,i,j,k,l,m,q,limit;
unsigned int flag,base,temp,tflag,sumdif,ps,pc,wrap;
unsigned int S[2],T[2],U[2],V[2],W[2],X[3],Y[4],Z[4];
int zflag;
unsigned int n=0;
FILE *Outfp;
Outfp = fopen("out6.dat","w");
/*********************************/
/* 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]);
tmpsav=tmpsiz;
limit=(tmptab[tmpsiz-1])>>16;
limit=limit*limit;
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
terror[0]=0; // clear error array
terror[1]=0;
terror[2]=0;
ps=p*p;
pc=ps*p;
zflag=0;
q=0;
count=0;
wrap=0;
for (h=offset; h<insize; h++) {
if (wrap>8) {
printf("input count=%d \n",h+1);
wrap=0;
}
else
wrap=wrap+1;
zloop:
if (zflag<2) {
d=input[3*h];
e=input[3*h+1];
c=input[3*h+2];
if (qflag==0) {
if (((d+e)/p)*p!=(d+e))
goto askip;
}
if (qflag==1) {
if (((d-e)/p)*p!=(d-e))
goto askip;
}
if (qflag==2) {
if ((d/p)*p!=d)
goto askip;
}
if (qflag==3) {
if ((e/p)*p!=e)
goto askip;
}
sumdif=0;
}
else {
d=input[3*(h-1)+1];
e=input[3*h+1];
c=input[3*h+2];
if (qflag==0) {
if (((d-e)/p)*p!=(d-e))
goto askip;
}
if (qflag==1) {
if (((d+e)/p)*p!=(d+e))
goto askip;
}
if (qflag==2) {
if ((d/p)*p!=d)
goto askip;
}
if (qflag==3) {
if ((e/p)*p!=e)
goto askip;
}
sumdif=1;
}
// dummy(d,e,0);
if ((d/p)*p==d) {
if (((d/2)*2==d)&&((d/ps)*ps!=d)) {
printf("error: p^2 does not divide d \n");
goto zskip;
}
tflag=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) {
if (((e/2)*2==e)&&((e/ps)*ps!=e)) {
printf("error: p^2 does not divide e \n");
goto zskip;
}
tflag=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==0)&&((((d+e)/2)*2==(d+e))&&(((d+e)/pc)*pc!=(d+e)))) {
printf("error: p^3 does not divide d+e \n");
goto zskip;
}
if ((sumdif==1)&&((((d+e)/2)*2==(d+e))&&(((d+e)/ps)*ps!=(d+e)))) {
printf("error: p^2 does not divide d-e \n");
goto zskip;
}
if (sumdif==0)
tflag=3;
else
tflag=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==0)&&((((d-e)/2)*2==(d-e))&&(((d-e)/ps)*ps!=(d-e)))) {
printf("error: p^2 does not divide d-e \n");
goto zskip;
}
if ((sumdif==1)&&((((d-e)/2)*2==(d-e))&&(((d-e)/pc)*pc!=(d-e)))) {
printf("error: p^3 does not divide d+e \n");
goto zskip;
}
if (sumdif==0)
tflag=2;
else
tflag=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;
}
}
if (select==0)
base=d;
if (select==1)
base=e;
if (select==2) {
if (sumdif==0)
base=d-e;
else
base=d+e;
}
if (select==3) {
if (sumdif==0)
base=d+e;
else
base=d-e;
}
if (pflag==1)
base=base*p;
if (pflag==2)
base=base*p*p;
if (select==4)
base=p;
if (select==5)
base=2;
if (select==6)
base=2*p;
/************************************/
/* 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==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);
}
S[0]=Y[2];
S[1]=Y[3];
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 (psflag==1) {
if (((k-1)/ps)*ps!=(k-1))
goto askip;
}
if (psflag==0) {
if (((k-1)/ps)*ps==(k-1))
goto askip;
}
bigresx(0, (k-1)/p, 0, k, U, base);
if ((U[0]!=0)||(U[1]!=1)) {
if (split==0) {
if (tflag==select)
terror[0]=8;
}
if (split==1) {
if ((tflag==0)&&(select==0))
terror[0]=8;
if ((tflag==1)&&(select==1))
terror[0]=8;
}
if (split==0) {
if ((select==5)&&((tflag==2)||(tflag==3)))
terror[0]=9;
}
goto askip;
}
aloop: S[0]=V[0];
S[1]=V[1];
save[m]=k;
if (m < savsiz) m=m+1;
else {
terror[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])) goto aloop;
}
/***********************************************/
/* output prime factors satisfying criterion */
/***********************************************/
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 (psflag==1) {
if (((i-1)/ps)*ps!=(i-1))
goto askip;
}
if (psflag==0) {
if (((i-1)/ps)*ps==(i-1))
goto askip;
}
bigresx(0, (i-1)/p, 0, i, U, base);
if ((U[0]!=0)||(U[1]!=1)) {
if (split==0) {
if (tflag==select)
terror[0]=8;
}
if (split==1) {
if ((tflag==0)&&(select==0))
terror[0]=8;
if ((tflag==1)&&(select==1))
terror[0]=8;
}
if (split==0) {
if ((select==5)&&((tflag==2)||(tflag==3)))
terror[0]=9;
}
goto askip;
}
if (T[0]<=limit) { // largest prime in table is 5801977
S[0]=T[0];
S[1]=T[1];
save[m]=i;
if (m < savsiz) m=m+1;
else {
terror[0]=3;
goto bskip;
}
goto cskip;
}
else {
terror[0]=4;
goto bskip;
}
}
}
}
cskip: if (psflag==1) {
T[0]=0;
T[1]=1;
differ(S,T);
quotient(T,T,ps);
bigprod(T[0],T[1],ps,X);
T[0]=0;
T[1]=1;
differ(S,T);
if ((X[1]!=T[0])||(X[2]!=T[1]))
goto askip;
}
if (psflag==0) {
T[0]=0;
T[1]=1;
differ(S,T);
quotient(T,T,ps);
bigprod(T[0],T[1],ps,X);
T[0]=0;
T[1]=1;
differ(S,T);
if ((X[1]==T[0])&&(X[2]==T[1]))
goto askip;
}
T[0]=0;
T[1]=1;
differ(S, T);
quotient(T, T, p);
bigresx(T[0], T[1], S[0], S[1], T, base);
if ((T[0]==0)&&(T[1]==1)) {
if (n+2>outsiz) {
terror[0]=6;
goto bskip;
}
output[n]=d;
output[n+1]=e;
output[n+2]=((m+1)<<16)|tflag;
if (m>1) {
if (q<999) {
terror[2*q+2]=d;
terror[2*q+3]=e;
}
q+=1;
terror[1]=q;
}
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];
}
if ((T[0]!=W[0]) || (T[1]!=W[1])) {
terror[0]=7;
goto bskip;
}
if ((prmflg==0)||(m==0)) {
n=n+3;
count=count+1;
}
}
else {
if (split==0) {
if (tflag==select)
terror[0]=8;
}
if (split==1) {
if ((tflag==0)&&(select==0))
terror[0]=8;
if ((tflag==1)&&(select==1))
terror[0]=8;
}
if (split==0) {
if ((select==5)&&((tflag==2)||(tflag==3)))
terror[0]=9;
}
goto askip;
}
}
else {
if (n+2>outsiz) {
terror[0]=6;
goto bskip;
}
output[n]=d;
output[n+1]=e;
output[n+2]=(m<<16)|tflag;
if (m>2) {
if (q<999) {
terror[2*q+2]=d;
terror[2*q+3]=e;
}
q+=1;
terror[1]=q;
}
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];
}
if ((S[0]!=W[0]) || (S[1]!=W[1])) {
terror[0]=7;
goto bskip;
}
if ((prmflg==0)||(m==1)) {
n=n+3;
count=count+1;
}
}
askip: if (zflag==2)
zflag=-1;
zflag+=1;
if (zflag==2)
goto zloop;
}
bskip:
output[n]=0xffffffff;
fprintf(Outfp," error0=%d \n",terror[0]);
fprintf(Outfp," count=%d \n",(n+1)/3);
for (i=0; i<(n+1)/3; i++)
fprintf(Outfp," %#10x %#10x %#10x \n",output[3*i],output[3*i+1],
output[3*i+2]);
zskip:
fclose(Outfp);
return(0);
}