/*****************************************************************************/
/* */
/* FACTOR (a**p-b**p)/(a-b) (when [(a**p+b**p)/(a+b)] is a pth power) */
/* 12/22/17 (dkc) */
/* */
/* Whether linear combinations of a and b are pth powers w.r.t. */
/* (a**p-b**p)/(a-b) is determined. */
/* */
/*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table9.h"
extern char *malloc();
unsigned int primed(unsigned int *out, unsigned int tsize,
unsigned int *table,unsigned int limit);
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 long long f);
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 differ(unsigned int *a, unsigned int *b);
unsigned int euclid(unsigned int d, unsigned int e);
int main ()
{
unsigned int p=3; // input prime (hardcoded to 3)
unsigned int numfac=2; // number of distinct prime factors of (a**p-b**p)/(a-b)
unsigned int split; // if set to 0, don't allow 2 and p to "split"
// if set to 1, only allow "split" 2 and p
// otherwise, allow both
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 offset=12500; // offset into the input look-up table (must be
// even and less than "insize"). Used to reduce
// execution time. insize=12810 for table9, (a,b) less than 5 million,
// and 7518 for table8, (a,b) between 5 and 10 million
unsigned int offset2=12500; // offset into the input look-up table (must be less than "insize"),
// for (f1,f2) coefficients, set to "offset" to check coverage
unsigned int mask2p;
//
unsigned int twopskip=1; // skip mask2p=4 if set (2p is a pth power w.r.t. (a^p-b^p)/(a-b), 2 and p not split)
unsigned int po2skip=1; // skip mask2p=8 if set (p/2 is a pth power w.r.t. (a^p-b^p)/(a-b), 2 and p not split)
unsigned int twoskip=1; // skip mask2p=2 if set (2 is a pth power w.r.t. (a^p-b^p)/(a-b), 2 and p not split)
unsigned int pskip=1; // skip mask2p=1 if set (p is a pth power w.r.t. (a^p-b^p)/(a-b), 2 and p not split)
unsigned int splitskip=1; // if set, skip mask2p=4 (2p is a pth power w.r.t. (a^p-b^p)/(a-b), split 2 and p)
unsigned int pa2skip=0; // skip mask2p=3 if set, for generators,
// 2 and p are pth powers w.r.t (a^p-b^p)/(a-b), 2 and p not split)
unsigned int nflag=0; // if set, use negative c4 value (normally set to 0)
unsigned int nocheck=0; // if set, coverage is not checked (normally set to 0)
//
unsigned int maskout=1; // output mask (normally set to 1)
unsigned int jump=1; // take jump if combination fails (normally set to 1)
unsigned int base2;
int c3;
int c4;
unsigned int pflag,pcount,qcount,zcount;
long long s3,s4;
unsigned int fail;
extern unsigned int input[];
extern unsigned int table[];
unsigned int *tmptab;
unsigned char *hits;
extern unsigned int output[];
extern unsigned int error[];
extern unsigned int count;
unsigned int tsize=1228;
unsigned int t2size=90000000;
unsigned int tmpsiz;
unsigned int outsiz=11000*7;
unsigned int save[20]; // solutions array
unsigned int savsiz=19; // size of solutions array minus one
unsigned int d,e,c,f;
unsigned int g,h,i,j,k,l,m,iters,xsave;
unsigned int flag,temp,sumdif,ps,twop;
unsigned int R[2],S[2],T[2],U[2],V[2],W[2],X[3],Y[4],Z[4];
int yflag,zflag;
unsigned int rflag,sflag,tflag,uflag,vflag,wflag,xflag;
unsigned int n;
FILE *Outfp;
Outfp = fopen("out32d.dat","w");
printf("largest (a,b) value= %d %d \n",input[3*offset],input[3*offset+1]);
if (offset>insize) {
printf("offset must be less than or equal to insize \n");
return(0);
}
if (offset!=(offset/2)*2) {
printf("offset must be even \n");
return(0);
}
if (offset2>insize) {
printf("offset2 must be less than or equal to insize \n");
return(0);
}
/*********************************/
/* extend prime look-up table */
/*********************************/
tmptab=(unsigned int*) malloc((t2size+1)*4);
if (tmptab==NULL) {
printf("not enough memory \n");
return(0);
}
tmpsiz=primed(tmptab, tsize, table, t2size);
j=0;
for (i=0; i<tmpsiz; i++) {
k=tmptab[i];
if ((k-1)==((k-1)/p)*p) {
tmptab[j]=k;
j=j+1;
}
}
tmpsiz=j;
printf("prime look-up table size=%d, last prime=%010x, %d \n",tmpsiz,tmptab[tmpsiz-1],tmptab[tmpsiz-1]);
/************************/
/* initialize hit table */
/************************/
hits=(unsigned char*) malloc(80000);
if (hits==NULL) {
printf("not enough memory \n");
return(0);
}
for (i=0; i<80000; i++)
hits[i]=0;
/**************/
/* mask loop */
/**************/
ps=p*p;
twop=p*2;
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
for (g=0; g<6; g++) {
if ((twopskip==0)&&(g==0)) {
mask2p=4;
base2=0;
split=0;
goto gskip;
}
if ((po2skip==0)&&(g==1)) {
mask2p=8;
base2=0;
split=0;
goto gskip;
}
if ((twoskip==0)&&(g==2)) {
mask2p=2;
base2=0;
split=0;
goto gskip;
}
if ((pskip==0)&&(g==3)) {
mask2p=1;
base2=1;
split=0;
goto gskip;
}
if ((splitskip==0)&&(g==4)) {
mask2p=4;
base2=0;
split=1;
goto gskip;
}
if ((pa2skip==0)&&(g==5)) {
mask2p=3;
split=0;
goto gskip;
}
continue;
/***********************************/
/* factor (d**p + e**p)/(d + e) */
/***********************************/
gskip:
count=0;
for (f=offset2; f<insize; f++) {
c3=(int)input[3*f];
c4=(int)input[3*f+1];
//printf("c3=%d c4=%d \n",c3,c4);
n=0;
zflag=0;
fail=0;
qcount=0;
zcount=0;
for (h=offset; h<insize; h++) {
zloop:
if (zflag<2) {
if (nflag!=0)
goto askip;
d=input[3*h];
e=input[3*h+1];
c=input[3*h+2];
sumdif=0;
}
else {
if (nflag==0)
goto askip;
d=input[3*(h-1)+1];
e=input[3*h+1];
c=input[3*h+2];
sumdif=1;
}
if ((d/p)*p==d) {
if (split==0) {
if ((d/2)*2!=d)
goto askip;
}
if (split==1) {
if ((d/2)*2==d)
goto askip;
}
yflag=0;
}
if ((e/p)*p==e) {
if (split==0) {
if ((e/2)*2!=e)
goto askip;
}
if (split==1) {
if ((e/2)*2==e)
goto askip;
}
yflag=1;
}
if (((d+e)/p)*p==(d+e)) {
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 (sumdif==0)
yflag=3;
else
yflag=2;
}
if (((d-e)/p)*p==(d-e)) {
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 (sumdif==0)
yflag=2;
else
yflag=3;
}
/************************************/
/* 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;
iters=0;
rflag=0xffffffff;
sflag=0xffffffff;
tflag=0xffffffff;
uflag=0xffffffff;
vflag=0xffffffff;
wflag=0xffffffff;
xflag=0xffffffff;
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;
}
if (iters<numfac) {
rflag=0;
if (zflag!=2) {
s3=(long long)c3*(long long)d+(long long)c4*(long long)e;
s4=(long long)c4*(long long)d+(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
else {
s3=(long long)c3*(long long)d-(long long)c4*(long long)e;
s4=(long long)c4*(long long)d-(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+512;
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+256;
if (g!=5) {
if (base2==0)
hugeres(0, (k-1)/p, 0, k, U, p*(unsigned long long)s3);
else
hugeres(0, (k-1)/p, 0, k, U, 2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+32;
if (base2==0)
hugeres(0, (k-1)/p, 0, k, U, p*(unsigned long long)s4);
else
hugeres(0, (k-1)/p, 0, k, U, 2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16;
if (base2==0)
hugeres(0, (k-1)/p, 0, k, U, p*p*(unsigned long long)s3);
else
hugeres(0, (k-1)/p, 0, k, U, 2*2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
if (base2==0)
hugeres(0, (k-1)/p, 0, k, U, p*p*(unsigned long long)s4);
else
hugeres(0, (k-1)/p, 0, k, U, 2*2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)p);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+4096;
R[0]=U[0];
R[1]=U[1];
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)2);
if ((U[0]==0)&&(U[1]==1)) {
rflag=rflag+8192;
if (split==1) {
if (((k-1)/ps)*ps!=(k-1))
error[0]=8;
}
}
if (g==5) {
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)d*(unsigned long long)d+(unsigned long long)e*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)(2*d)*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
if (g!=5) {
if ((R[0]==U[0])&&(R[1]==U[1]))
rflag=rflag+32768;
hugeres(0, (k-1)/p, 0, k, U, (unsigned long long)(2*p));
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16384;
}
}
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])) goto aloop;
iters=iters+1;
if (iters<numfac) {
xflag=wflag;
wflag=vflag;
vflag=uflag;
uflag=tflag;
tflag=sflag;
sflag=rflag;
}
else {
if ((S[0]!=0)||(S[1]!=1))
goto askip;
}
}
/***********************************************/
/* 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+=twop) {
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;
}
if (iters<numfac) {
rflag=0;
if (zflag!=2) {
s3=(long long)c3*(long long)d+(long long)c4*(long long)e;
s4=(long long)c4*(long long)d+(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
else {
s3=(long long)c3*(long long)d-(long long)c4*(long long)e;
s4=(long long)c4*(long long)d-(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+512;
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+256;
if (g!=5) {
if (base2==0)
hugeres(0, (i-1)/p, 0, i, U, p*(unsigned long long)s3);
else
hugeres(0, (i-1)/p, 0, i, U, 2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+32;
if (base2==0)
hugeres(0, (i-1)/p, 0, i, U, p*(unsigned long long)s4);
else
hugeres(0, (i-1)/p, 0, i, U, 2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16;
if (base2==0)
hugeres(0, (i-1)/p, 0, i, U, p*p*(unsigned long long)s3);
else
hugeres(0, (i-1)/p, 0, i, U, 2*2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
if (base2==0)
hugeres(0, (i-1)/p, 0, i, U, p*p*(unsigned long long)s4);
else
hugeres(0, (i-1)/p, 0, i, U, 2*2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)p);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+4096;
R[0]=U[0];
R[1]=U[1];
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)2);
if ((U[0]==0)&&(U[1]==1)) {
rflag=rflag+8192;
if (split==1) {
if (((i-1)/ps)*ps!=(i-1))
error[0]=8;
}
}
if (g==5) {
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)d*(unsigned long long)d+(unsigned long long)e*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)(2*d)*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
if (g!=5) {
if ((R[0]==U[0])&&(R[1]==U[1]))
rflag=rflag+32768;
hugeres(0, (i-1)/p, 0, i, U, (unsigned long long)(2*p));
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16384;
}
}
iters=iters+1;
if (iters<numfac) {
xflag=wflag;
wflag=vflag;
vflag=uflag;
uflag=tflag;
tflag=sflag;
sflag=rflag;
}
if (T[0]<=4) { // largest prime in table is 0x270eb
S[0]=T[0];
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: if ((S[0]==0)&&(S[1]==1)) {
if (iters!=numfac)
goto askip;
else
goto dskip;
}
if ((S[0]==0)&&(S[1]==save[m])) {
if (iters==numfac)
goto dskip;
else
goto askip;
}
if (iters!=(numfac-1))
goto askip;
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);
rflag=0;
if (zflag!=2) {
s3=(long long)c3*(long long)d+(long long)c4*(long long)e;
s4=(long long)c4*(long long)d+(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
else {
s3=(long long)c3*(long long)d-(long long)c4*(long long)e;
s4=(long long)c4*(long long)d-(long long)c3*(long long)e;
if (s3<0)
s3=-s3;
if (s4<0)
s4=-s4;
}
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+512;
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+256;
if (g!=5) {
if (base2==0)
hugeres(T[0], T[1], S[0], S[1], U, p*(unsigned long long)s3);
else
hugeres(T[0], T[1], S[0], S[1], U, 2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+32;
if (base2==0)
hugeres(T[0], T[1], S[0], S[1], U, p*(unsigned long long)s4);
else
hugeres(T[0], T[1], S[0], S[1], U, 2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16;
if (base2==0)
hugeres(T[0], T[1], S[0], S[1], U, p*p*(unsigned long long)s3);
else
hugeres(T[0], T[1], S[0], S[1], U, 2*2*(unsigned long long)s3);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
if (base2==0)
hugeres(T[0], T[1], S[0], S[1], U, p*p*(unsigned long long)s4);
else
hugeres(T[0], T[1], S[0], S[1], U, 2*2*(unsigned long long)s4);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)p);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+4096;
R[0]=U[0];
R[1]=U[1];
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)2);
if ((U[0]==0)&&(U[1]==1)) {
rflag=rflag+8192;
if (split==1) {
V[0]=0;
V[1]=1;
differ(S,V);
quotient(V,V,ps);
bigprod(V[0],V[1],ps,X);
V[0]=0;
V[1]=1;
differ(S,V);
if ((X[1]!=V[0])||(X[2]!=V[1]))
error[0]=8;
}
}
if (g==5) {
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)d*(unsigned long long)d+(unsigned long long)e*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+2;
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)(2*d)*(unsigned long long)e);
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+1;
}
if (g!=5) {
if ((R[0]==U[0])&&(R[1]==U[1]))
rflag=rflag+32768;
hugeres(T[0], T[1], S[0], S[1], U, (unsigned long long)(2*p));
if ((U[0]==0)&&(U[1]==1))
rflag=rflag+16384;
}
dskip: if (n+6>outsiz) {
error[0]=6;
goto bskip;
}
xflag=rflag&sflag&tflag&uflag&vflag&wflag&xflag;
if ((xflag>>12)==mask2p) {
pcount=0;
if (g!=5)
pflag=xflag&0x333;
else {
xsave=xflag&3;
pflag=xflag&0x330;
}
while (pflag!=0) {
if ((pflag&1)==1)
pcount=pcount+1;
pflag=pflag>>1;
}
if ((pcount!=0)&&(pcount!=2)) {
fail=fail+1;
if (jump!=0)
goto eskip;
}
if (pcount!=0) {
if ((g==5)&&(xsave!=3)) {
printf("error: d=%d e=%d xsave=%d \n",d,e,xsave);
fail=fail+1;
goto eskip;
}
qcount=qcount+1;
}
else
zcount=zcount+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])) {
error[0]=7;
goto bskip;
}
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]))
error[1]+=1;
count=count+1;
if ((numfac==2)&&(split==0)&&(psflag==0)) {
if ((yflag==0)||(yflag==1)) {
if (rflag!=sflag)
error[0]=8;
T[0]=0;
T[1]=1;
differ(W,T);
quotient(T,T,ps);
bigprod(T[0],T[1],ps,X);
T[0]=0;
T[1]=1;
differ(W,T);
if ((X[1]!=T[0])||(X[2]!=T[1]))
error[0]=8;
}
}
}
//
// S[0]=0, S[1]=1
//
else {
if (iters!=numfac)
goto askip;
if (n+6>outsiz) {
error[0]=6;
goto bskip;
}
xflag=rflag&sflag&tflag&uflag&vflag&wflag&xflag;
if ((xflag>>12)==mask2p) {
pcount=0;
if (g!=5)
pflag=xflag&0x333;
else {
xsave=xflag&3;
pflag=xflag&0x330;
}
while (pflag!=0) {
if ((pflag&1)==1)
pcount=pcount+1;
pflag=pflag>>1;
}
if ((pcount!=0)&&(pcount!=2)) {
fail=fail+1;
if (jump!=0)
goto eskip;
}
if (pcount!=0) {
if ((g==5)&&(xsave!=3)) {
printf("error: d=%d e=%d xsave=%d \n",d,e,xsave);
fail=fail+1;
goto eskip;
}
qcount=qcount+1;
}
else
zcount=zcount+1;
}
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])) {
error[0]=7;
goto bskip;
}
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]))
error[1]+=1;
count=count+1;
if ((numfac==2)&&(split==0)&&(psflag==0)) {
if ((yflag==0)||(yflag==1)) {
if (rflag!=sflag)
error[0]=8;
T[0]=0;
T[1]=1;
differ(W,T);
quotient(T,T,ps);
bigprod(T[0],T[1],ps,X);
T[0]=0;
T[1]=1;
differ(W,T);
if ((X[1]!=T[0])||(X[2]!=T[1]))
error[0]=8;
}
}
}
askip: if (zflag==2)
zflag=-1;
zflag+=1;
if (zflag==2)
goto zloop;
} // h loop
if (fail!=0)
printf("fail=%d, f=%d, g=%d \n",fail,c3,c4);
if ((qcount!=0)&&(fail==0)) {
printf("passed: f=%d, g=%d, count=%d, zcount=%d \n",c3,c4,qcount,zcount);
fprintf(Outfp,"passed: f=%d, g=%d, count=%d, zcount=%d \n",c3,c4,qcount,zcount);
output[n]=c3;
output[n+1]=c4;
output[n+2]=qcount;
n=n+3;
hits[f]=hits[f]+1;
}
eskip:
i=0;
} // f loop
} // g loop
if (nocheck==0) {
printf("checking coverage \n");
qcount=0;
for (f=offset2; f<insize; f++) {
qcount=qcount+1;
if (nflag==0) {
if (hits[f]==0)
printf("index=%d, miss=%d %d \n",f,input[3*f],input[3*f+1]);
}
else {
if (hits[f]==0)
printf("index=%d, miss=%d %d \n",f,input[3*f],input[3*f+4]);
}
}
printf("count=%d \n",qcount);
}
bskip:
printf("error=%d \n",error[0]);
fclose(Outfp);
return(0);
}
/******************************************************************************
* *
* 128-BIT DIFFERENCE *
* 11/13/06 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void bigbigd(unsigned int *a, unsigned int *b) {
unsigned int a0,a1,a2,a3,b0,b1,b2,b3;
unsigned int s0,s1,s2,s3,temp,c,c1,c2,c3;
a0=*a;
a1=*(a+1);
a2=*(a+2);
a3=*(a+3);
b0=*b;
b1=*(b+1);
b2=*(b+2);
b3=*(b+3);
b0=~b0;
b1=~b1;
b2=~b2;
b3=~b3;
temp=b3+1;
c=carry(b3,1,temp);
b3=temp;
temp=b2+c;
c=carry(b2,c,temp);
b2=temp;
temp=b1+c;
c=carry(b1,c,temp);
b1=temp;
b0=b0+c;
s3=a3+b3;
c3=carry(a3,b3,s3);
s2=a2+b2;
c2=carry(a2,b2,s2);
s1=a1+b1;
c1=carry(a1,b1,s1);
s0=a0+b0;
temp=s2+c3;
c2+=carry(s2,c3,temp);
s2=temp;
temp=s1+c2;
c1+=carry(s1,c2,temp);
s1=temp;
s0=s0+c1;
*b=s0;
*(b+1)=s1;
*(b+2)=s2;
*(b+3)=s3;
return;
}
/*****************************************************************************/
/* */
/* 128/64 DIVISION (UNSIGNED) */
/* 11/14/06 (dkc) */
/* */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
unsigned int lmbd(unsigned int mode, unsigned int a);
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 i,d0,d1,dshift,ashift,count,flag;
unsigned int shift,c,c0,c1,c2,temp,temp0,temp1,temp2,temp3;
if (d2==0) {
if ((a0==0)&&(a1==0)&&(a2==0)&&(a3<d3)) {
*quotient=0;
*(quotient+1)=0;
*(quotient+2)=0;
*(quotient+3)=0;
return;
}
}
else {
if ((a0==0)&&(a1==0)&&(a2<d2)) {
*quotient=0;
*(quotient+1)=0;
*(quotient+2)=0;
*(quotient+3)=0;
return;
}
}
dshift=lmbd(1,d2);
if (d2==0)
dshift+=lmbd(1,d3);
dshift+=64;
ashift=lmbd(1,a0);
if (a0==0)
ashift+=lmbd(1,a1);
if ((a0|a1)==0)
ashift+=lmbd(1,a2);
if ((a0|a1|a2)==0)
ashift+=lmbd(1,a3);
shift=dshift-ashift;
count=shift+1;
d0=0;
d1=0;
if (shift<32) {
if (shift!=0) {
d1=d2>>(32-shift);
d2=(d2<<shift)|(d3>>(32-shift));
}
else
d1=0; // added to get MSVC to work
d3=d3<<shift;
flag=3;
shift=32-shift;
}
else {
shift=shift-32;
d1=d2;
d2=d3;
d3=0;
if (shift<32) {
if (shift!=0) {
d0=d1>>(32-shift);
d1=(d1<<shift)|(d2>>(32-shift));
}
else
d0=0; // added to get MSVC to work
d2=d2<<shift;
flag=2;
shift=32-shift;
}
else {
shift=shift-32;
d0=d1;
d1=d2;
d2=0;
if (shift<32) {
if (shift!=0) // added to get MSVC to work
d0=(d0<<shift)|(d1>>(32-shift));
d1=d1<<shift;
flag=1;
shift=32-shift;
}
else {
shift=shift-32;
d0=d1;
d1=0;
d0=d0<<shift;
flag=0;
shift=32-shift;
}
}
}
d0=~d0;
d1=~d1;
d2=~d2;
d3=~d3;
temp=d3+1;
c=carry(d3,1,temp);
d3=temp;
temp=d2+c;
c=carry(d2,c,temp);
d2=temp;
temp=d1+c;
c=carry(d1,c,temp);
d1=temp;
d0=d0+c;
for (i=0; i<count; i++) {
temp3=a3+d3;
c2=carry(a3,d3,temp3);
temp=a2+c2;
c1=carry(a2,c2,temp);
temp2=temp+d2;
c1+=carry(temp,d2,temp2);
temp=a1+c1;
c0=carry(a1,c1,temp);
temp1=temp+d1;
c0+=carry(temp,d1,temp1);
temp0=a0+d0+c0;
if ((temp0>>31)==0) {
a0=temp0<<1;
if ((temp1>>31)!=0)
c=1;
else
c=0;
a0=a0+c;
a1=temp1<<1;
if ((temp2>>31)!=0)
c=1;
else
c=0;
a1=a1+c;
a2=temp2<<1;
if ((temp3>>31)!=0)
c=1;
else
c=0;
a2=a2+c;
a3=temp3<<1;
a3=a3+1;
}
else {
a0=a0<<1;
if ((a1>>31)!=0)
c=1;
else
c=0;
a0=a0+c;
a1=a1<<1;
if ((a2>>31)!=0)
c=1;
else
c=0;
a1=a1+c;
a2=a2<<1;
if ((a3>>31)!=0)
c=1;
else
c=0;
a2=a2+c;
a3=a3<<1;
}
}
shift=shift-1;
if (flag==3) {
a0=0;
a1=0;
a2=0;
a3=a3<<shift;
a3=a3>>shift;
}
else {
if (flag==2) {
a0=0;
a1=0;
a2=a2<<shift;
a2=a2>>shift;
}
else {
if (flag==1) {
a0=0;
a1=a1<<shift;
a1=a1>>shift;
}
else {
a0=a0<<shift;
a0=a0>>shift;
}
}
}
*quotient=a0;
*(quotient+1)=a1;
*(quotient+2)=a2;
*(quotient+3)=a3;
return;
}
/******************************************************************************
* *
* 128-BIT SUM *
* 11/13/06 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void bigbigs(unsigned int *a, unsigned int *b) {
unsigned int a0,a1,a2,a3,b0,b1,b2,b3;
unsigned int s0,s1,s2,s3,temp,c1,c2,c3;
a0=*a;
a1=*(a+1);
a2=*(a+2);
a3=*(a+3);
b0=*b;
b1=*(b+1);
b2=*(b+2);
b3=*(b+3);
s3=a3+b3;
c3=carry(a3,b3,s3);
s2=a2+b2;
c2=carry(a2,b2,s2);
s1=a1+b1;
c1=carry(a1,b1,s1);
s0=a0+b0;
temp=s2+c3;
c2+=carry(s2,c3,temp);
s2=temp;
temp=s1+c2;
c1+=carry(s1,c2,temp);
s1=temp;
s0=s0+c1;
*b=s0;
*(b+1)=s1;
*(b+2)=s2;
*(b+3)=s3;
return;
}
/******************************************************************************
* *
* 64x64 MULTIPLY (UNSIGNED) *
* 11/14/06 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void bigbigp(unsigned int a0, unsigned int a2, unsigned int m0,
unsigned int m2, unsigned int *product) {
unsigned int a1,a3,m1,m3,temp;
unsigned int p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15;
unsigned int s0,s1,s2,s3;
unsigned int c1,c2,c3;
a1=a0&0xffff;
a0=a0>>16;
a3=a2&0xffff;
a2=a2>>16;
m1=m0&0xffff;
m0=m0>>16;
m3=m2&0xffff;
m2=m2>>16;
p0=a0*m0;
p1=a0*m1;
p2=a1*m0;
p3=a0*m2;
p4=a1*m1;
p5=a2*m0;
p6=a0*m3;
p7=a1*m2;
p8=a2*m1;
p9=a3*m0;
p10=a1*m3;
p11=a2*m2;
p12=a3*m1;
p13=a2*m3;
p14=a3*m2;
p15=a3*m3;
s3=p15+(p14<<16);
c3=carry(p15,(p14<<16),s3);
temp=s3+(p13<<16);
c3+=carry(s3,(p13<<16),temp);
s3=temp;
s2=p12+(p14>>16);
c2=carry(p12,(p14>>16),s2);
temp=s2+(p13>>16);
c2+=carry(s2,(p13>>16),temp);
s2=temp;
temp=s2+p11;
c2+=carry(s2,p11,temp);
s2=temp;
temp=s2+p10;
c2+=carry(s2,p10,temp);
s2=temp;
temp=s2+(p9<<16);
c2+=carry(s2,(p9<<16),temp);
s2=temp;
temp=s2+(p8<<16);
c2+=carry(s2,(p8<<16),temp);
s2=temp;
temp=s2+(p7<<16);
c2+=carry(s2,(p7<<16),temp);
s2=temp;
temp=s2+(p6<<16);
c2+=carry(s2,(p6<<16),temp);
s2=temp;
s1=p5+(p9>>16);
c1=carry(p5,(p9>16),s1);
temp=s1+(p8>>16);
c1+=carry(s1,(p8>>16),temp);
s1=temp;
temp=s1+(p7>>16);
c1+=carry(s1,(p7>>16),temp);
s1=temp;
temp=s1+(p6>>16);
c1+=carry(s1,(p6>>16),temp);
s1=temp;
temp=s1+p4;
c1+=carry(s1,p4,temp);
s1=temp;
temp=s1+p3;
c1+=carry(s1,p3,temp);
s1=temp;
temp=s1+(p2<<16);
c1+=carry(s1,(p2<<16),temp);
s1=temp;
temp=s1+(p1<<16);
c1+=carry(s1,(p1<<16),temp);
s1=temp;
s0=p0+(p2>>16);
s0=s0+(p1>>16);
temp=s2+c3;
c2+=carry(s2,c3,temp);
s2=temp;
temp=s1+c2;
c1+=carry(s1,c2,temp);
s1=temp;
s0=s0+c1;
*product=s0;
*(product+1)=s1;
*(product+2)=s2;
*(product+3)=s3;
return;
}
/******************************************************************************
* *
* 64x64 MULTIPLY (UNSIGNED) *
* 03/23/12 (dkc) *
* *
******************************************************************************/
void bignewp(unsigned int a0, unsigned int a1, unsigned int b0,
unsigned int b1, unsigned int *product) {
unsigned long long ac,bc,ad,bd,mid34,upper,lower;
bd=(unsigned long long)a1*(unsigned long long)b1;
ad=(unsigned long long)a0*(unsigned long long)b1;
bc=(unsigned long long)b0*(unsigned long long)a1;
ac=(unsigned long long)a0*(unsigned long long)b0;
mid34=(bd>>32)+(bc&0xffffffff)+(ad&0xffffffff);
upper=ac+(bc>>32)+(ad>>32)+(mid34>>32);
lower=(mid34<<32)|(bd&0xffffffff);
*product=(unsigned int)(upper>>32);
*(product+1)=(unsigned int)(upper&0xffffffff);
*(product+2)=(unsigned int)(lower>>32);
*(product+3)=(unsigned int)(lower&0xffffffff);
return;
}
/******************************************************************************
* *
* 64x32 MULTIPLY (UNSIGNED) *
* 11/14/06 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void bigprod(unsigned int a0, unsigned int a2, unsigned int m0,
unsigned int *product) {
unsigned int a1,a3,m1,temp;
unsigned int p0,p1,p2,p3,p4,p5,p6,p7;
unsigned int s1,s2,s3;
unsigned int c2,c3;
a1=a0&0xffff;
a0=a0>>16;
a3=a2&0xffff;
a2=a2>>16;
m1=m0&0xffff;
m0=m0>>16;
p0=a0*m0;
p1=a0*m1;
p2=a1*m0;
p3=a1*m1;
p4=a2*m0;
p5=a2*m1;
p6=a3*m0;
p7=a3*m1;
s3=p7+(p6<<16);
c3=carry(p7,(p6<<16),s3);
temp=s3+(p5<<16);
c3+=carry(s3,(p5<<16),temp);
s3=temp;
s2=p4+(p6>>16);
c2=carry(p4,(p6>>16),s2);
temp=s2+(p5>>16);
c2+=carry(s2,(p5>>16),temp);
s2=temp;
temp=s2+p3;
c2+=carry(s2,p3,temp);
s2=temp;
temp=s2+(p2<<16);
c2+=carry(s2,(p2<<16),temp);
s2=temp;
temp=s2+(p1<<16);
c2+=carry(s2,(p1<<16),temp);
s2=temp;
s1=p0+(p2>>16);
s1=s1+(p1>>16);
temp=s2+c3;
c2+=carry(s2,c3,temp);
s2=temp;
s1=s1+c2;
*product=s1;
*(product+1)=s2;
*(product+2)=s3;
return;
}
/******************************************************************************/
/* */
/* FIND LEAST RESIDUE (UNSIGNED) */
/* 11/13/06 (dkc) */
/* */
/* This C subroutine finds the least residue of "base**exponent" modulus q. */
/* The binary representation of the exponent is used to efficiently compute */
/* the least residue. */
/* */
/******************************************************************************/
void differ(unsigned int *a, unsigned int *b);
void bigbigp(unsigned int d0, unsigned int d1, unsigned int m0, unsigned int m1,
unsigned int *output);
void bigbigq(unsigned int a0, unsigned int a1, unsigned int a2, unsigned int a3,
unsigned int *quotient, unsigned int d0, unsigned int d1);
void bigbigd(unsigned int *subtrahend, unsigned int *minuend);
void quotient(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor);
void bigprod(unsigned int a0, unsigned int a1, unsigned int m0,
unsigned int *output);
void bigresx(unsigned int exp0, unsigned int exp1, unsigned int q0,
unsigned int q1, unsigned int *output, unsigned int base) {
int i,index;
unsigned int save[128];
unsigned int B[2],E[2],F[2],T[2],P[4],R[4],X[3];
B[0]=0; // load base
B[1]=base; // load base
if ((q0==0)&&(base>q1)) { // check for base greater than q
quotient(B, T, q1); // base/q1
bigprod(T[0], T[1], q1, X); // (base/q1)*q1
T[0]=X[1];
T[1]=X[2];
differ(B, T); // base=base-(base/q1)*q1
B[0]=T[0];
B[1]=T[1];
}
if (base==1) // check for base equal to 1
goto askip;
E[0]=0; // load exponent
E[1]=1; // load exponent
F[0]=exp0; // load input exponent
F[1]=exp1; // load input exponent
for (i=0; i<64; i++) {
save[2*i]=B[0]; // save base
save[2*i+1]=B[1]; // save base
bigprod(E[0], E[1], 2, X); // exponent=exponent*2
if (X[1]>F[0]) { // compare exponent to input exponent
index=i;
break;
}
if ((X[1]==F[0])&&(X[2]>F[1])) { // compare exponent to input exponent
index=i;
break;
}
E[0]=X[1];
E[1]=X[2];
bigbigp(B[0], B[1], B[0], B[1], P); // base=base**2
bigbigq(P[0], P[1], P[2], P[3], R, q0, q1); // base/q
bigbigp(R[2], R[3], q0, q1, R); // (base/q)*q
bigbigd(P, R); // base=base-(base/q)*q
B[0]=R[2];
B[1]=R[3];
if ((E[0]==F[0])&&(E[1]==F[1])) // compare exponent to input exponent
goto askip;
}
B[0]=save[2*index]; // load base
B[1]=save[2*index+1]; // load base
T[0]=E[0];
T[1]=E[1];
differ(F, T); // input exponent-=exponent
F[0]=T[0];
F[1]=T[1];
for (i=index-1; i>=0; i--) {
quotient(E, E, 2); // exponent=exponent/2
T[0]=E[0];
T[1]=E[1];
differ(F, T); // input exponent-=exponent
if ((T[0]&0x80000000)==0) { // check if non-negative
F[0]=T[0];
F[1]=T[1];
bigbigp(B[0], B[1], save[2*i], save[2*i+1], P); // base*=save[i]
bigbigq(P[0], P[1], P[2], P[3], R, q0, q1); // base/q
bigbigp(R[2], R[3], q0, q1, R); // (base/q)*q
bigbigd(P, R); // base=base-(base/p)*q
B[0]=R[2];
B[1]=R[3];
}
if ((F[0]==0)&&(F[1]==0)) // check if input exponent is zero
break;
}
askip:
*output=B[0]; // store least residue
*(output+1)=B[1]; // store least residue
}
/*****************************************************************************/
/* */
/* CARRY */
/* 11/15/06 (dkc) */
/* */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum) {
unsigned int c;
if ((a&b)>>31!=0)
c=1;
else {
if ((a|b)>>31==0)
c=0;
else {
if (sum>=0x80000000)
c=0;
else
c=1;
}
}
return c;
}
/******************************************************************************
* *
* COMPUTE 64-BIT DIFFERENCE *
* 11/13/06 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void differ(unsigned int *a, unsigned int *b) {
unsigned int high0,high1,low0,low1,templo,temphi,temp,c;
high0=*a;
low0=*(a+1);
high1=*b;
low1=*(b+1);
low1=~low1;
high1=~high1;
temp=low1+1;
c=carry(low1,1,temp);
low1=temp;
high1+=c;
templo=low0+low1;
c=carry(low0,low1,templo);
temphi=high0+high1+c;
*b=temphi;
*(b+1)=templo;
return;
}
/*****************************************************************************/
/* */
/* EUCLIDEAN G.C.D. */
/* 11/14/06 (dkc) */
/* */
/*****************************************************************************/
unsigned int euclid(unsigned int d, unsigned int e) {
unsigned int a,b,temp;
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;
return(a);
}
/******************************************************************************/
/* */
/* FIND LEAST RESIDUE (UNSIGNED) */
/* 11/13/06 (dkc) */
/* */
/* This C subroutine finds the least residue of "base**exponent" modulus q. */
/* The binary representation of the exponent is used to efficiently compute */
/* the least residue. */
/* */
/******************************************************************************/
void differ(unsigned int *a, unsigned int *b);
void bignewp(unsigned int d0, unsigned int d1, unsigned int m0, unsigned int m1,
unsigned int *output);
void bigbigq(unsigned int a0, unsigned int a1, unsigned int a2, unsigned int a3,
unsigned int *quotient, unsigned int d0, unsigned int d1);
void bigbigd(unsigned int *subtrahend, unsigned int *minuend);
void bigprod(unsigned int a0, unsigned int a1, unsigned int m0,
unsigned int *output);
void quotient(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor);
void hugeres(unsigned int exp0, unsigned int exp1, unsigned int q0,
unsigned int q1, unsigned int *output, unsigned long long base) {
int i,index;
unsigned int save[128];
unsigned int B[2],E[2],F[2],T[2],P[4],R[4],X[3];
unsigned int base0,base1;
base0=(unsigned int)(base>>32);
base1=(unsigned int)(base&0xffffffff);
B[0]=base0; // load base
B[1]=base1; // load base
if ((base0>q0)||((base0==q0)&&(base1>q1))) { // check for base greater than q
bigbigq(0, 0, B[0], B[1], R, q0, q1); // base/q
bignewp(R[2], R[3], q0, q1, R); // (base/q)*q
T[0]=R[2];
T[1]=R[3];
differ(B, T); // base=base-(base/q)*q
B[0]=T[0];
B[1]=T[1];
}
if ((B[0]==0)&&(B[1]==1)) // check for base equal to 1
goto askip;
E[0]=0; // load exponent
E[1]=1; // load exponent
F[0]=exp0; // load input exponent
F[1]=exp1; // load input exponent
for (i=0; i<64; i++) {
save[2*i]=B[0]; // save base
save[2*i+1]=B[1]; // save base
bigprod(E[0], E[1], 2, X); // exponent=exponent*2
if (X[1]>F[0]) { // compare exponent to input exponent
index=i;
break;
}
if ((X[1]==F[0])&&(X[2]>F[1])) { // compare exponent to input exponent
index=i;
break;
}
E[0]=X[1];
E[1]=X[2];
bignewp(B[0], B[1], B[0], B[1], P); // base=base**2
bigbigq(P[0], P[1], P[2], P[3], R, q0, q1); // base/q
bignewp(R[2], R[3], q0, q1, R); // (base/q)*q
bigbigd(P, R); // base=base-(base/q)*q
B[0]=R[2];
B[1]=R[3];
if ((E[0]==F[0])&&(E[1]==F[1])) // compare exponent to input exponent
goto askip;
}
B[0]=save[2*index]; // load base
B[1]=save[2*index+1]; // load base
T[0]=E[0];
T[1]=E[1];
differ(F, T); // input exponent-=exponent
F[0]=T[0];
F[1]=T[1];
for (i=index-1; i>=0; i--) {
quotient(E, E, 2); // exponent=exponent/2
T[0]=E[0];
T[1]=E[1];
differ(F, T); // input exponent-=exponent
if ((T[0]&0x80000000)==0) { // check if non-negative
F[0]=T[0];
F[1]=T[1];
bignewp(B[0], B[1], save[2*i], save[2*i+1], P); // base*=save[i]
bigbigq(P[0], P[1], P[2], P[3], R, q0, q1); // base/q
bignewp(R[2], R[3], q0, q1, R); // (base/q)*q
bigbigd(P, R); // base=base-(base/p)*q
B[0]=R[2];
B[1]=R[3];
}
if ((F[0]==0)&&(F[1]==0)) // check if input exponent is zero
break;
}
askip:
*output=B[0]; // store least residue
*(output+1)=B[1]; // store least residue
}
/*****************************************************************************/
/* */
/* LEFT-MOST BIT DETECTION */
/* 11/15/06 (dkc) */
/* */
/*****************************************************************************/
unsigned int lmbd(unsigned int mode, unsigned int a) {
unsigned int i,mask,count;
if (mode==0)
a=~a;
if (a==0)
return(32);
count=0;
mask=0x80000000;
for (i=0; i<32; i++) {
if ((a&mask)!=0)
break;
count+=1;
mask>>=1;
}
return(count);
}
/*****************************************************************************/
/* */
/* COMPUTE PRIMES */
/* 11/04/15 (dkc) */
/* */
/*****************************************************************************/
#include <math.h>
unsigned int primed(unsigned int *out, unsigned int tsize,
unsigned int *table,unsigned int limit) {
unsigned int d;
unsigned int i,j,k,l,flag,count;
count=tsize;
for (i=0; i<tsize; i++)
out[i]=table[i];
j=table[tsize-1]+1;
/********************/
/* loop through d */
/********************/
for (d=j; d<=limit; 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;
/************************************************/
/* look for prime factors using look-up table */
/************************************************/
l=(unsigned int)(10.0+sqrt((double)d));
k=0;
if (l>table[tsize-1])
return(0);
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)
out[count]=d;
count=count+flag;
}
return(count);
}
/*****************************************************************************/
/* */
/* 64/32 UNSIGNED DIVIDE */
/* 11/14/06 (dkc) */
/* */
/*****************************************************************************/
unsigned int lmbd(unsigned int mode, unsigned int a);
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void quotient(unsigned int *dividend, unsigned int *quotient,
unsigned int divisor) {
unsigned int i;
unsigned int shift,dshift,ashift,d0,d1,div0,div1,c,temp0,temp1,t0,t1;
d0=*dividend;
d1=*(dividend+1);
if ((d0==0)&&(d1<divisor)) {
*quotient=0;
*(quotient+1)=0;
return;
}
dshift=lmbd(1, divisor);
dshift=dshift+32;
ashift=lmbd(1, d0);
if (d0==0)
ashift+=lmbd(1, d1);
shift=dshift-ashift;
if (shift<32) {
div1=divisor<<shift;
if (shift!=0)
div0=divisor>>(32-shift);
else // added to get MSVC to work
div0=0; //
}
else {
if (shift!=32)
div0=divisor<<(shift-32);
else // added to get MSVC to work
div0=divisor;
div1=0;
}
t0=~div0;
t1=~div1;
temp1=t1+1;
c=carry(t1,1,temp1);
t1=temp1;
t0=t0+c;
shift+=1;
for (i=0; i<shift; i++) {
temp1=d1+t1;
c=carry(d1,t1,temp1);
temp0=d0+t0+c;
if ((temp0>>31)==0) {
d0=temp0<<1;
if ((temp1>>31)!=0)
c=1;
else
c=0;
d0=d0+c;
d1=temp1<<1;
d1=d1+1;
}
else {
d0=d0<<1;
if ((d1>>31)!=0)
c=1;
else
c=0;
d0=d0+c;
d1=d1<<1;
}
}
if (shift>32) {
d0=d0<<(64-shift);
d0=d0>>(64-shift);
}
else {
d0=0;
if (shift!=32) { // added to get MSVC to work
d1=d1<<(32-shift);
d1=d1>>(32-shift);
}
}
*quotient=d0;
*(quotient+1)=d1;
return;
}