﻿ big least residue
/******************************************************************************/
/*									      */
/*  FIND LEAST RESIDUE (UNSIGNED)					      */
/*  01/24/14 (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];

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
for (i=0; i<64; i++) {
save[2*i]=B[0];		      // save base
save[2*i+1]=B[1];		      // save base
X[0]=(E[0]<<1)|(E[1]>>31);	      // exponent=exponent*2
X[1]=E[1]<<1;
if (F[0]==0) {
if ((F[1]&0x80000000)==0) {
if (X[1]>F[1]) {   // compare exponent to input exponent
index=i;
break;
}
}
else {
if (X[1]==0) {
index=i;
break;
}
}
}
else {
if (X[0]>F[0]) {
index=i;
break;
}
}
E[0]=X[0];
E[1]=X[1];
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
}
F[1]=F[1]^E[1];
for (i=index-1; i>=0; i--) {
E[1]=(E[1]>>1)|(E[0]<<31);	      // exponent=exponent/2
E[0]=E[0]>>1;
T[0]=F[0]&E[0];		      // check for bit
T[1]=F[1]&E[1];
if ((T[0]!=0)||(T[1]!=0)) {	      // check if 0
F[1]=F[1]^E[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/q)*q
B[0]=R[2];
B[1]=R[3];
}
if ((F[0]==0)&&(F[1]==0))	      // check if input exponent is zero
break;
}