﻿ huge residue
```/******************************************************************************/
/*									      */
/*  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];

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
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
}
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;
}