/******************************************************************************/
/*									      */
/*  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];

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
   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
      goto askip;
   }
B[0]=save[2*index];		      // load base
B[1]=save[2*index+1];		      // load base
F[0]=F[0]^E[0]; 		      // mask bit
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[0]=F[0]^E[0];		      // mask bit
      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;
   }
askip:
*output=B[0];			      // store least residue
*(output+1)=B[1];		      // store least residue
}