/******************************************************************************/
/*									      */
/*  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 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 int base0,
	     unsigned int base1) {
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]=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
   bigbigp(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];
   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
}