/*****************************************************************************/ /* */ /* FACTOR (a**p+b**p)/(a+b) */ /* 11/07/06 (dkc) */ /* */ /* This C program determines if (a**p+b**p)/(a+b) is prime or p times a */ /* prime. p must divide a, b, a-b, or a+b if p>3. p must divide a, b, */ /* a-b or p**2 must divide a+b if p=3. If p divides a, then q must divide */ /* a. If p divides b, then q must divide b. If p divides a-b or a+b, */ /* then q must divide a-b or a+b. Whether q is a pth power modulus */ /* [(a**p+b**p)/(a+b)] is determined. p should not divide q. */ /* */ /* The output is "(a<<16)|b". If q is not a pth power modulus [(a**p+b**p)/*/ /* (a+b)], then an error is indicated ("error[1]" is set to a non-zero */ /* value). */ /* */ /*****************************************************************************/ #include <math.h> #include <stdio.h> #include "table9.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 bigresx(unsigned int a, unsigned int b, unsigned int c, unsigned int d, unsigned int *e, unsigned int f); 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 5000. // The maximum "dbeg" value for p=5 is about 1000. // The maximum "dbeg" value for p=7 is about 250. // The maximum "dbeg" value for p=11 is about 50. // unsigned int p=5; // input prime unsigned int dbeg=500; // starting "a" value unsigned int dend=1; // ending "a" value //unsigned int stop=1; 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=51; // q value 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 tsize=172; unsigned int tmpsiz; unsigned int outsiz=1999; unsigned int d,e,a,b,temp; unsigned int i,j,k,l,m; unsigned int flag; unsigned int S[2],T[2],U[2],X[3],Y[4],Z[4]; unsigned int ps,pb; unsigned int n=0; double sqrt2=1.4142135; FILE *Outfp; Outfp = fopen("out19.dat","w"); /*********************************/ /* extend prime look-up table */ /*********************************/ for (i=0; i<tsize; i++) tmptab[i] = (int)(table[i]); tmpsiz=tsize; for (d=1033; d<170000; 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; /************************************************/ /* look for prime factors using 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; 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) */ /***********************************/ pb=p*base; ps=p*p; error[0]=0; // clear error array error[1]=0; error[2]=0; error[3]=0; error[4]=0; error[5]=0; count=0; for (d=dbeg; d>=dend; d--) { for (e=d-1; e>0; e--) { // if (e!=stop) continue; /******************************************/ /* check if p divides d, e, d+e or d-e */ /******************************************/ if ((d/pb)*pb==d) goto zskip; if ((e/pb)*pb==e) goto zskip; if (p==3) { if (sumdif==1) { if (((d+e)/p)*p==(d+e)) { if (((d+e)/ps)*ps!=(d+e)) continue; } } else { if (((d-e)/p)*p==(d-e)) { if (((d-e)/ps)*ps!=(d-e)) continue; } } } if (((d+e)/p)*p==(d+e)) { if (((d+e)/base)*base==(d+e)) goto zskip; if (((d-e)/base)*base==(d-e)) goto zskip; else continue; } if (((d-e)/p)*p==(d-e)) { if (((d-e)/base)*base==(d-e)) goto zskip; if (((d+e)/base)*base==(d+e)) goto zskip; else continue; } else 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; /***********************/ /* 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) */ /************************************/ 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); 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); 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]; /************************************************/ /* look for prime factors using look-up table */ /************************************************/ if (S[0]==0) l = 32 - lmbd(1, S[1]); else l = 64 - lmbd(1, S[0]); j=l-(l/2)*2; l=l/2; 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; } } for (i=0; i<=k; i++) { m = tmptab[i]; quotient(S, T, m); bigprod(T[0], T[1], m, X); if ((S[0]==X[1]) && (S[1]==X[2])) goto askip; } /***********************************************/ /* output prime factors satisfying criterion */ /***********************************************/ if (flag==1) { for (i=tmptab[tmpsiz-1]; i<=l; i+=2) { quotient(S, T, i); bigprod(T[0], T[1], i, X); if ((X[1]==S[0]) && (X[2]==S[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], U, base); if ((U[0]!=0)||(U[1]!=1)) { error[1]+=1; goto askip; } if (n>outsiz) { error[0]=6; n=n-1; } output[n]=((int)(d) << 16) | (int)(e); n=n+1; count=count+1; askip:dummy(d,e,8); } } bskip: output[n]=-1; fprintf(Outfp," error0=%d error1=%d \n", error[0],error[1]); fprintf(Outfp," count=%d \n",n); if (n!=0) { for (i=0; i<n-1; i++) fprintf(Outfp," %#10x \n",output[i]); } if (error[1]!=0) printf(" error \n"); fclose(Outfp); return(0); }