--- pgp/src/mpilib.c	2018/04/24 16:40:00	1.1.1.4
+++ pgp/src/mpilib.c	2018/04/24 16:42:05	1.1.1.6
@@ -7,13 +7,14 @@
 	Boulder, CO 80304
 	(303) 541-0140
 
-	(c) Copyright 1986-92 by Philip Zimmermann.  All rights reserved.
+	(c) Copyright 1986-1994 by Philip Zimmermann.  All rights reserved.
 	The author assumes no liability for damages resulting from the use
 	of this software, even if the damage results from defects in this
 	software.  No warranty is expressed or implied.  The use of this
 	cryptographic software for developing weapon systems is expressly
 	forbidden.
 
+
 	These routines implement all of the multiprecision arithmetic
 	necessary for number-theoretic cryptographic algorithms such as
 	ElGamal, Diffie-Hellman, Rabin, or factoring studies for large
@@ -139,25 +140,25 @@ boolean mp_addc
 	(register unitptr r1,register unitptr r2,register boolean carry)
 	/* multiprecision add with carry r2 to r1, result in r1 */
 	/* carry is incoming carry flag-- value should be 0 or 1 */
-{	register unit x;
+{
+	register unit x;
 	short precision;	/* number of units to add */
 	precision = global_precision;
 	make_lsbptr(r1,precision);
 	make_lsbptr(r2,precision);
-	while (precision--)
-	{
-		if (carry)
-		{	x = *r1 + *r2 + 1;
+	while (precision--) {
+		if (carry) {
+			x = *r1 + *r2 + 1;
 			carry = (*r2  >= (unit)(~ *r1));
-		} else
-		{	x = *r1 + *r2;
+		} else {
+			x = *r1 + *r2;
 			carry = (x < *r1) ;
 		}
 		post_higherunit(r2);
 		*post_higherunit(r1) = x;
 	}
-	return(carry);		/* return the final carry flag bit */
-}	/* mp_addc */
+	return carry;		/* return the final carry flag bit */
+} /* mp_addc */
 #endif  /* mp_addc */
 
 #ifndef mp_subb
@@ -165,44 +166,45 @@ boolean mp_subb
 	(register unitptr r1,register unitptr r2,register boolean borrow)
 	/* multiprecision subtract with borrow, r2 from r1, result in r1 */
 	/* borrow is incoming borrow flag-- value should be 0 or 1 */
-{	register unit x;
+{
+	register unit x;
 	short precision;	/* number of units to subtract */
 	precision = global_precision;
 	make_lsbptr(r1,precision);
 	make_lsbptr(r2,precision);
-	while (precision--)
-	{	if (borrow)
-		{	x = *r1 - *r2 - 1;
+	while (precision--) {
+		if (borrow) {
+			x = *r1 - *r2 - 1;
 			borrow = (*r1 <= *r2);
-		} else
-		{	x = *r1 - *r2;
+		} else {
+			x = *r1 - *r2;
 			borrow = (*r1 < *r2);
 		}
 		post_higherunit(r2);
 		*post_higherunit(r1) = x;
 	}
-	return(borrow);	/* return the final carry/borrow flag bit */
-}	/* mp_subb */
+	return borrow;	/* return the final carry/borrow flag bit */
+} /* mp_subb */
 #endif  /* mp_subb */
 
 #ifndef mp_rotate_left
 boolean mp_rotate_left(register unitptr r1,register boolean carry)
 	/* multiprecision rotate left 1 bit with carry, result in r1. */
 	/* carry is incoming carry flag-- value should be 0 or 1 */
-{	register int precision;	/* number of units to rotate */
+{
+	register int precision;	/* number of units to rotate */
 	unsigned int mcarry = carry, nextcarry; /* int is supposed to be
 						 * the efficient size for ops*/
 	precision = global_precision;
 	make_lsbptr(r1,precision);
-	while (precision--)
-	{
+	while (precision--) {
 		nextcarry = (((signedunit) *r1) < 0);
 		*r1 = (*r1 << 1) | mcarry;
 		mcarry = nextcarry;
 		pre_higherunit(r1);
 	}
-	return(nextcarry);	/* return the final carry flag bit */
-}	/* mp_rotate_left */
+	return nextcarry;	/* return the final carry flag bit */
+} /* mp_rotate_left */
 #endif  /* mp_rotate_left */
 
 /************** end of primitives that should be in assembly *************/
@@ -216,7 +218,8 @@ void mp_shift_right_bits(register unitpt
 	/*	multiprecision shift right bits, result in r1.
 		bits is how many bits to shift, must be <= UNITSIZE.
 	*/
-{	register short precision;	/* number of units to shift */
+{
+	register short precision;	/* number of units to shift */
 	register unit carry,nextcarry,bitmask;
 	register short unbits;
 	if (bits==0) return;	/* shift zero bits is a no-op */
@@ -233,84 +236,90 @@ void mp_shift_right_bits(register unitpt
 			pre_lowerunit(r1);
 		}
 	} else {
-		while (precision--)
-		{	nextcarry = *r1 & bitmask;
+		while (precision--) {
+			nextcarry = *r1 & bitmask;
 			*r1 >>= bits;
 			*r1 |= carry << unbits;
 			carry = nextcarry;
 			pre_lowerunit(r1);
 		}
 	}
-}	/* mp_shift_right_bits */
+} /* mp_shift_right_bits */
 
 
 #ifndef mp_compare
 short mp_compare(register unitptr r1,register unitptr r2)
-/*	Compares multiprecision integers *r1, *r2, and returns:
-	-1 iff *r1 < *r2
-	 0 iff *r1 == *r2
-	+1 iff *r1 > *r2
-*/
-{	register short precision;	/* number of units to compare */
+/*
+ * Compares multiprecision integers *r1, *r2, and returns:
+ * -1 iff *r1 < *r2
+ *  0 iff *r1 == *r2
+ * +1 iff *r1 > *r2
+ */
+{
+	register short precision;	/* number of units to compare */
 
 	precision = global_precision;
 	make_msbptr(r1,precision);
 	make_msbptr(r2,precision);
-	do
-	{	if (*r1 < *r2)
-			return(-1);
+	do {
+		if (*r1 < *r2)
+			return -1;
 		if (*post_lowerunit(r1) > *post_lowerunit(r2))
-			return(1);
+			return 1;
 	} while (--precision);
-	return(0);	/*  *r1 == *r2  */
-}	/* mp_compare */
+	return 0;	/*  *r1 == *r2  */
+} /* mp_compare */
 #endif /* mp_compare */
 
 
 boolean mp_inc(register unitptr r)
-	/* Increment multiprecision integer r. */
-{	register short precision;
+/* Increment multiprecision integer r. */
+{
+	register short precision;
 	precision = global_precision;
 	make_lsbptr(r,precision);
-	do
-	{	if ( ++(*r) ) return(0);	/* no carry */
+	do {
+		if ( ++(*r) ) return 0;	/* no carry */
 		post_higherunit(r);
 	} while (--precision);
-	return(1);		/* return carry set */
-}	/* mp_inc */
+	return 1;		/* return carry set */
+} /* mp_inc */
 
 
 boolean mp_dec(register unitptr r)
-	/* Decrement multiprecision integer r. */
-{	register short precision;
+/* Decrement multiprecision integer r. */
+{
+	register short precision;
 	precision = global_precision;
 	make_lsbptr(r,precision);
-	do
-	{	if ( (signedunit) (--(*r)) != (signedunit) -1 )
-			return(0);	/* no borrow */
+	do {
+		if ( (signedunit) (--(*r)) != (signedunit) -1 )
+			return 0;	/* no borrow */
 		post_higherunit(r);
 	} while (--precision);
-	return(1);		/* return borrow set */
-}	/* mp_dec */
+	return 1;		/* return borrow set */
+} /* mp_dec */
 
 
 void mp_neg(register unitptr r)
 	/* Compute 2's complement, the arithmetic negative, of r */
-{	register short precision;	/* number of units to negate */
+{
+	register short precision;	/* number of units to negate */
 	precision = global_precision;
 	mp_dec(r);	/* 2's complement is 1's complement plus 1 */
-	do	/* now do 1's complement */
-	{	*r = ~(*r);
+	do {	/* now do 1's complement */
+		*r = ~(*r);
 		r++;
 	} while (--precision);
-}	/* mp_neg */
+} /* mp_neg */
 
 #ifndef mp_move
 void mp_move(register unitptr dst,register unitptr src)
-{	register short precision;	/* number of units to move */
+{
+	register short precision;	/* number of units to move */
 	precision = global_precision;
 	do { *dst++ = *src++; } while (--precision);
-}	/* mp_move */
+} /* mp_move */
 #endif /* mp_move */
 
 void mp_init(register unitptr r, word16 value)
@@ -323,38 +332,42 @@ void mp_init(register unitptr r, word16
 #ifdef UNIT8
 	*post_higherunit(r) = value >> UNITSIZE;
 #endif
-}	/* mp_init */
+} /* mp_init */
 
 
 short significance(register unitptr r)
 	/* Returns number of significant units in r */
-{	register short precision;
+{
+	register short precision;
 	precision = global_precision;
 	make_msbptr(r,precision);
-	do
-	{	if (*post_lowerunit(r))
-			return(precision);
+	do {
+		if (*post_lowerunit(r))
+			return precision;
 	} while (--precision);
-	return(precision);
-}	/* significance */
+	return precision;
+} /* significance */
 
 
 #ifndef unitfill0
 void unitfill0(unitptr r,word16 unitcount)
 	/* Zero-fill the unit buffer r. */
-{	while (unitcount--) *r++ = 0;
-}	/* unitfill0 */
+{
+	while (unitcount--) *r++ = 0;
+} /* unitfill0 */
 #endif  /* unitfill0 */
 
 
 int mp_udiv(register unitptr remainder,register unitptr quotient,
 	register unitptr dividend,register unitptr divisor)
 	/* Unsigned divide, treats both operands as positive. */
-{	int bits;
+{
+	int bits;
 	short dprec;
 	register unit bitmask;
+
 	if (testeq(divisor,0))
-		return(-1);	/* zero divisor means divide error */
+		return -1;	/* zero divisor means divide error */
 	mp_init0(remainder);
 	mp_init0(quotient);
 	/* normalize and compute number of bits in dividend first */
@@ -363,15 +376,15 @@ int mp_udiv(register unitptr remainder,r
 	rescale(quotient,global_precision,dprec);
 	make_msbptr(quotient,dprec);
 
-	while (bits--)
-	{	mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
-		if (mp_compare(remainder,divisor) >= 0)
-		{	mp_sub(remainder,divisor);
+	while (bits--) {
+		mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
+		if (mp_compare(remainder,divisor) >= 0) {
+			mp_sub(remainder,divisor);
 			stuff_bit(quotient,bitmask);
 		}
 		bump_2bitsniffers(dividend,quotient,bitmask);
 	}
-	return(0);
+	return 0;
 } /* mp_udiv */
 
 
@@ -380,12 +393,13 @@ int mp_udiv(register unitptr remainder,r
 
 int mp_recip(register unitptr quotient,register unitptr divisor)
 	/* Compute reciprocal (quotient) as 1/divisor.  Used by faster modmult. */
-{	int bits;
+{
+	int bits;
 	short qprec;
 	register unit bitmask;
 	unit remainder[MAX_UNIT_PRECISION];
 	if (testeq(divisor,0))
-		return(-1);	/* zero divisor means divide error */
+		return -1;	/* zero divisor means divide error */
 	mp_init0(remainder);
 	mp_init0(quotient);
 
@@ -398,16 +412,16 @@ int mp_recip(register unitptr quotient,r
 	rescale(quotient,global_precision,qprec);
 	make_msbptr(quotient,qprec);
 
-	while (bits--)
-	{	mp_shift_left(remainder);
-		if (mp_compare(remainder,divisor) >= 0)
-		{	mp_sub(remainder,divisor);
+	while (bits--) {
+		mp_shift_left(remainder);
+		if (mp_compare(remainder,divisor) >= 0) {
+			mp_sub(remainder,divisor);
 			stuff_bit(quotient,bitmask);
 		}
 		bump_bitsniffer(quotient,bitmask);
 	}
 	mp_init0(remainder);	/* burn sensitive data left on stack */
-	return(0);
+	return 0;
 } /* mp_recip */
 #endif	/* UPTON_OR_SMITH */
 
@@ -415,7 +429,8 @@ int mp_recip(register unitptr quotient,r
 int mp_div(register unitptr remainder,register unitptr quotient,
 	register unitptr dividend,register unitptr divisor)
 	/* Signed divide, either or both operands may be negative. */
-{	boolean dvdsign,dsign;
+{
+	boolean dvdsign,dsign;
 	int status;
 	dvdsign = (boolean)(mp_tstminus(dividend)!=0);
 	dsign = (boolean)(mp_tstminus(divisor)!=0);
@@ -424,26 +439,28 @@ int mp_div(register unitptr remainder,re
 	status = mp_udiv(remainder,quotient,dividend,divisor);
 	if (dvdsign) mp_neg(dividend);		/* repair caller's dividend */
 	if (dsign) mp_neg(divisor);		/* repair caller's divisor */
-	if (status<0) return(status);		/* divide error? */
+	if (status<0) return status;		/* divide error? */
 	if (dvdsign) mp_neg(remainder);
 	if (dvdsign ^ dsign) mp_neg(quotient);
-	return(status);
+	return status;
 } /* mp_div */
 
 
 word16 mp_shortdiv(register unitptr quotient,
 	register unitptr dividend,register word16 divisor)
-/*	This function does a fast divide and mod on a multiprecision dividend
-	using a short integer divisor returning a short integer remainder.
-	This is an unsigned divide.  It treats both operands as positive.
-	It is used mainly for faster printing of large numbers in base 10.
-*/
-{	int bits;
+/*
+ * This function does a fast divide and mod on a multiprecision dividend
+ * using a short integer divisor returning a short integer remainder.
+ * This is an unsigned divide.  It treats both operands as positive.
+ * It is used mainly for faster printing of large numbers in base 10.
+ */
+{
+	int bits;
 	short dprec;
 	register unit bitmask;
 	register word16 remainder;
 	if (!divisor)	/* if divisor == 0 */
-		return(-1);	/* zero divisor means divide error */
+		return -1;	/* zero divisor means divide error */
 	remainder=0;
 	mp_init0(quotient);
 	/* normalize and compute number of bits in dividend first */
@@ -452,65 +469,67 @@ word16 mp_shortdiv(register unitptr quot
 	rescale(quotient,global_precision,dprec);
 	make_msbptr(quotient,dprec);
 
-	while (bits--)
-	{	remainder <<= 1;
+	while (bits--) {
+		remainder <<= 1;
 		if (sniff_bit(dividend,bitmask))
 			remainder++;
-		if (remainder >= divisor)
-		{	remainder -= divisor;
+		if (remainder >= divisor) {
+			remainder -= divisor;
 			stuff_bit(quotient,bitmask);
 		}
 		bump_2bitsniffers(dividend,quotient,bitmask);
 	}
-	return(remainder);
+	return remainder;
 } /* mp_shortdiv */
 
 
 int mp_mod(register unitptr remainder,
 	register unitptr dividend,register unitptr divisor)
-	/* Unsigned divide, treats both operands as positive. */
-{	int bits;
+/* Unsigned divide, treats both operands as positive. */
+{
+	int bits;
 	short dprec;
 	register unit bitmask;
 	if (testeq(divisor,0))
-		return(-1);	/* zero divisor means divide error */
+		return -1;	/* zero divisor means divide error */
 	mp_init0(remainder);
 	/* normalize and compute number of bits in dividend first */
 	init_bitsniffer(dividend,bitmask,dprec,bits);
 
-	while (bits--)
-	{	mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
+	while (bits--) {
+		mp_rotate_left(remainder,(boolean)(sniff_bit(dividend,bitmask)!=0));
 		msub(remainder,divisor);
 		bump_bitsniffer(dividend,bitmask);
 	}
-	return(0);
+	return 0;
 } /* mp_mod */
 
 
 word16 mp_shortmod(register unitptr dividend,register word16 divisor)
-/*	This function does a fast mod operation on a multiprecision dividend
-	using a short integer modulus returning a short integer remainder.
-	This is an unsigned divide.  It treats both operands as positive.
-	It is used mainly for fast sieve searches for large primes.
-*/
+/*
+ * This function does a fast mod operation on a multiprecision dividend
+ * using a short integer modulus returning a short integer remainder.
+ * This is an unsigned divide.  It treats both operands as positive.
+ * It is used mainly for fast sieve searches for large primes.
+ */
 {	int bits;
 	short dprec;
 	register unit bitmask;
 	register word16 remainder;
 	if (!divisor)	/* if divisor == 0 */
-		return(-1);	/* zero divisor means divide error */
+		return -1;	/* zero divisor means divide error */
 	remainder=0;
 	/* normalize and compute number of bits in dividend first */
 	init_bitsniffer(dividend,bitmask,dprec,bits);
 
-	while (bits--)
-	{	remainder <<= 1;
+	while (bits--) {
+		remainder <<= 1;
 		if (sniff_bit(dividend,bitmask))
 			remainder++;
 		if (remainder >= divisor) remainder -= divisor;
 		bump_bitsniffer(dividend,bitmask);
 	}
-	return(remainder);
+	return remainder;
 } /* mp_shortmod */
 
 
@@ -520,27 +539,30 @@ word16 mp_shortmod(register unitptr divi
 
 int mp_mult(register unitptr prod,
 	register unitptr multiplicand, register unitptr multiplier)
-	/*	Computes multiprecision prod = multiplicand * multiplier */
-{	/*	Uses interleaved comb multiply algorithm.
-		This improved multiply more than twice as fast as a Russian
-		peasant multiply, because it does a lot fewer shifts.
-		Must have global_precision set to the size of the multiplicand
-		plus UNITSIZE-1 SLOP_BITS.  Produces a product that is the sum
-		of the lengths of the multiplier and multiplicand.
-
-		BUG ALERT:  Unfortunately, this code has a bug.  It fails for
-		some numbers.  One such example:
-		x=   59DE 60CE 2345 8091 A02B 2A1C DBC3 8BE5
-		x*x= 59DE 60CE 2345 26B3 993B 67A5 2499 0B7D
-		     52C8 CDC7 AFB3 61C8 243C 741B
-		--which is obviously wrong.  The answer should be:
-		x*x= 1F8C 607B 5EA6 C061 2714 04A9 A0C6 A17A
-		     C9AB 6095 C62F 3756 3843 E4D0 3950 7AD9
-		We'll have to fix this some day.  In the meantime, we'll
-		just have the compiler skip over this code.
-
-		BUG NOTE: Possibly fixed.  Needs testing.
-	*/
+/*
+ * Computes multiprecision prod = multiplicand * multiplier
+ *
+ * Uses interleaved comb multiply algorithm.
+ * This improved multiply more than twice as fast as a Russian
+ * peasant multiply, because it does a lot fewer shifts.
+ * Must have global_precision set to the size of the multiplicand
+ * plus UNITSIZE-1 SLOP_BITS.  Produces a product that is the sum
+ * of the lengths of the multiplier and multiplicand.
+ *
+ * BUG ALERT:  Unfortunately, this code has a bug.  It fails for
+ * some numbers.  One such example:
+ * x=   59DE 60CE 2345 8091 A02B 2A1C DBC3 8BE5
+ * x*x= 59DE 60CE 2345 26B3 993B 67A5 2499 0B7D
+ *      52C8 CDC7 AFB3 61C8 243C 741B
+ * --which is obviously wrong.  The answer should be:
+ * x*x= 1F8C 607B 5EA6 C061 2714 04A9 A0C6 A17A
+ *      C9AB 6095 C62F 3756 3843 E4D0 3950 7AD9
+ * We'll have to fix this some day.  In the meantime, we'll
+ * just have the compiler skip over this code.
+ *
+ * BUG NOTE: Possibly fixed.  Needs testing.
+ */
+{
 	int bits;
 	register unit bitmask;
 	unitptr product, mplier, temp;
@@ -551,33 +573,33 @@ int mp_mult(register unitptr prod,
 	mp_init(prod+tohigher(global_precision),0);
 
 	if (testeq(multiplicand,0))
-		return(0);	/* zero multiplicand means zero product */
+		return 0;	/* zero multiplicand means zero product */
 
 	mp_move(mplicand,multiplicand);	/* save it from damage */
 
 	normalize(multiplier,mprec);
 	if (!mprec)
-		return(0);
+		return 0;
 
 	make_lsbptr(multiplier,mprec);
 	bitmask = 1;	/* start scan at LSB of multiplier */
 
-	do	/* UNITSIZE times */
-	{	/* do for bits 0-15 */
+	do {	/* UNITSIZE times */
+		/* do for bits 0-15 */
 		product = prod;
 		mplier = multiplier;
 		mprec2 = mprec;
-		while (mprec2--)	/* do for each word in multiplier */
-		{
-			if (sniff_bit(mplier,bitmask))
-			{	if (mp_addc(product,multiplicand,0))	/* ripple carry */
+		while (mprec2--) {	/* do for each word in multiplier */
+		
+			if (sniff_bit(mplier,bitmask)) {
+				if (mp_addc(product,multiplicand,0))	/* ripple carry */
 				{	/* After 1st time thru, this is rarely encountered. */
 					temp = msbptr(product,global_precision);
 					pre_higherunit(temp);
 					/* temp now points to LSU of carry region. */
 					unmake_lsbptr(temp,global_precision);
 					mp_inc(temp);
-				}	/* ripple carry */
+				} /* ripple carry */
 			}
 			pre_higherunit(mplier);
 			pre_higherunit(product);
@@ -590,8 +612,8 @@ int mp_mult(register unitptr prod,
 
 	mp_move(multiplicand,mplicand);	/* recover */
 
-	return(0);	/* normal return */
-}	/* mp_mult */
+	return 0;	/* normal return */
+} /* mp_mult */
 
 #endif	/* COMB_MULT */
 
@@ -603,50 +625,54 @@ int mp_mult(register unitptr prod,
 
 int mp_mult(register unitptr prod,
 	register unitptr multiplicand,register unitptr multiplier)
-	/* Computes multiprecision prod = multiplicand * multiplier */
-{	/* Uses "Russian peasant" multiply algorithm. */
+/*
+ * Computes multiprecision prod = multiplicand * multiplier
+ * Uses "Russian peasant" multiply algorithm.
+ */
+{
 	int bits;
 	register unit bitmask;
 	short mprec;
 	mp_init(prod,0);
 	if (testeq(multiplicand,0))
-		return(0);	/* zero multiplicand means zero product */
+		return 0;	/* zero multiplicand means zero product */
 	/* normalize and compute number of bits in multiplier first */
 	init_bitsniffer(multiplier,bitmask,mprec,bits);
 
-	while (bits--)
-	{	mp_shift_left(prod);
+	while (bits--) {
+		mp_shift_left(prod);
 		if (sniff_bit(multiplier,bitmask))
 			mp_add(prod,multiplicand);
 		bump_bitsniffer(multiplier,bitmask);
 	}
-	return(0);
-}	/* mp_mult */
+	return 0;
+} /* mp_mult */
 
 
 
-/*	mp_modmult computes a multiprecision multiply combined with a
-	modulo operation.  This is the most time-critical function in
-	this multiprecision arithmetic library for performing modulo
-	exponentiation.  We experimented with different versions of modmult,
-	depending on the machine architecture and performance requirements.
-	We will either use a Russian Peasant modmult (peasant_modmult), 
-	Charlie Merritt's modmult (merritt_modmult), Jimmy Upton's
-	modmult (upton_modmult), or Thad Smith's modmult (smith_modmult).
-	On machines with a hardware atomic multiply instruction,
-	Smith's modmult is fastest.  It can utilize assembly subroutines to
-	speed up the hardware multiply logic and trial quotient calculation.
-	If the machine lacks a fast hardware multiply, Merritt's modmult
-	is preferred, which doesn't call any assembly multiply routine.
-	We use the alias names mp_modmult, stage_modulus, and modmult_burn
-	for the corresponding true names, which depend on what flavor of
-	modmult we are using.
-
-	Before making the first call to mp_modmult, you must set up the
-	modulus-dependant precomputated tables by calling stage_modulus.
-	After making all the calls to mp_modmult, you call modmult_burn to
-	erase the tables created by stage_modulus that were left in memory.
-*/
+/*
+ * mp_modmult computes a multiprecision multiply combined with a
+ * modulo operation.  This is the most time-critical function in
+ * this multiprecision arithmetic library for performing modulo
+ * exponentiation.  We experimented with different versions of modmult,
+ * depending on the machine architecture and performance requirements.
+ * We will either use a Russian Peasant modmult (peasant_modmult), 
+ * Charlie Merritt's modmult (merritt_modmult), Jimmy Upton's
+ * modmult (upton_modmult), or Thad Smith's modmult (smith_modmult).
+ * On machines with a hardware atomic multiply instruction,
+ * Smith's modmult is fastest.  It can utilize assembly subroutines to
+ * speed up the hardware multiply logic and trial quotient calculation.
+ * If the machine lacks a fast hardware multiply, Merritt's modmult
+ * is preferred, which doesn't call any assembly multiply routine.
+ * We use the alias names mp_modmult, stage_modulus, and modmult_burn
+ * for the corresponding true names, which depend on what flavor of
+ * modmult we are using.
+ * 
+ * Before making the first call to mp_modmult, you must set up the
+ * modulus-dependant precomputated tables by calling stage_modulus.
+ * After making all the calls to mp_modmult, you call modmult_burn to
+ * erase the tables created by stage_modulus that were left in memory.
+ */
 
 #ifdef COUNTMULTS
 /* "number of modmults" counters, used for performance studies. */
@@ -661,56 +687,61 @@ static unsigned int tally_modsquares = 0
 static unitptr pmodulus = 0;	/* used only by mp_modmult */
 
 int stage_peasant_modulus(unitptr n)
-/*	Must pass modulus to stage_modulus before calling modmult.
-	Assumes that global_precision has already been adjusted to the
-	size of the modulus, plus SLOP_BITS.
-*/
+/*
+ * Must pass modulus to stage_modulus before calling modmult.
+ * Assumes that global_precision has already been adjusted to the
+ * size of the modulus, plus SLOP_BITS.
+ */
 {	/* For this simple version of modmult, just copy unit pointer. */
 	pmodulus = n;
-	return(0);	/* normal return */
-}	/* stage_peasant_modulus */
+	return 0;	/* normal return */
+} /* stage_peasant_modulus */
 
 
 int peasant_modmult(register unitptr prod,
 	unitptr multiplicand,register unitptr multiplier)
-{	/*	"Russian peasant" multiply algorithm, combined with a modulo
-		operation.  This is a simple naive replacement for Merritt's
-		faster modmult algorithm.  References global unitptr "modulus".
-		Computes:  prod = (multiplicand*multiplier) mod modulus
-		WARNING: All the arguments must be less than the modulus!
-	*/
+/* "Russian peasant" multiply algorithm, combined with a modulo
+ * operation.  This is a simple naive replacement for Merritt's
+ * faster modmult algorithm.  References global unitptr "modulus".
+ * Computes:  prod = (multiplicand*multiplier) mod modulus
+ * WARNING: All the arguments must be less than the modulus!
+ */
+{
 	int bits;
 	register unit bitmask;
 	short mprec;
 	mp_init(prod,0);
 /*	if (testeq(multiplicand,0))
-		return(0); */	/* zero multiplicand means zero product */
+		return 0; */	/* zero multiplicand means zero product */
 	/* normalize and compute number of bits in multiplier first */
 	init_bitsniffer(multiplier,bitmask,mprec,bits);
 
-	while (bits--)
-	{	mp_shift_left(prod);
+	while (bits--) {
+		mp_shift_left(prod);
 		msub(prod,pmodulus);	/* turns mult into modmult */
-		if (sniff_bit(multiplier,bitmask))
-		{	mp_add(prod,multiplicand);
+		if (sniff_bit(multiplier,bitmask)) {
+			mp_add(prod,multiplicand);
 			msub(prod,pmodulus);	/* turns mult into modmult */
 		}
 		bump_bitsniffer(multiplier,bitmask);
 	}
-	return(0);
-}	/* peasant_modmult */
+	return 0;
+} /* peasant_modmult */
 
 
-/*	If we are using a version of mp_modmult that uses internal tables
-	in memory, we have to call modmult_burn() at the end of mp_modexp.
-	This is so that no sensitive data is left in memory after the program
-	exits.  The Russian peasant method doesn't use any such tables.
-*/
+/*
+ * If we are using a version of mp_modmult that uses internal tables
+ * in memory, we have to call modmult_burn() at the end of mp_modexp.
+ * This is so that no sensitive data is left in memory after the program
+ * exits.  The Russian peasant method doesn't use any such tables.
+ */
 void peasant_burn(void)
-/*	Alias for modmult_burn, called only from mp_modexp().  Destroys
-	internal modmult tables.  This version does nothing because no
-	tables are used by the Russian peasant modmult. */
-{ }	/* peasant_burn */
+/*
+ * Alias for modmult_burn, called only from mp_modexp().  Destroys
+ * internal modmult tables.  This version does nothing because no
+ * tables are used by the Russian peasant modmult.
+ */
+{ } /* peasant_burn */
 
 #endif	/* PEASANT */
 
@@ -732,7 +763,8 @@ void peasant_burn(void)
 
 static void mp_lshift_unit(register unitptr r1)
 /*	Shift r1 1 whole unit to the left.  Used only by modmult function. */
-{	register short precision;
+{
+	register short precision;
 	register unitptr r2;
 	precision = global_precision;
 	make_msbptr(r1,precision);
@@ -740,68 +772,74 @@ static void mp_lshift_unit(register unit
 	while (--precision)
 		*post_lowerunit(r1) = *pre_lowerunit(r2);
 	*r1 = 0;
-}	/* mp_lshift_unit */
+} /* mp_lshift_unit */
 
 
 /* moduli_buf contains shifted images of the modulus, set by stage_modulus */
 static unit moduli_buf[UNITSIZE-1][MAX_UNIT_PRECISION] = {0};
 static unitptr moduli[UNITSIZE] = /* contains pointers into moduli_buf */
-{	0
-	,&moduli_buf[ 0][0], &moduli_buf[ 1][0], &moduli_buf[ 2][0], &moduli_buf[ 3][0],
-	 &moduli_buf[ 4][0], &moduli_buf[ 5][0], &moduli_buf[ 6][0]
+{	 0,              moduli_buf[ 0], moduli_buf[ 1], moduli_buf[ 2],
+	 moduli_buf[ 3], moduli_buf[ 4], moduli_buf[ 5], moduli_buf[ 6]
 #ifndef UNIT8
-	,&moduli_buf[ 7][0]
-	,&moduli_buf[ 8][0], &moduli_buf[ 9][0], &moduli_buf[10][0], &moduli_buf[11][0], 
-	 &moduli_buf[12][0], &moduli_buf[13][0], &moduli_buf[14][0]
+	,moduli_buf[ 7] ,moduli_buf[ 8], moduli_buf[ 9], moduli_buf[10],
+	 moduli_buf[11], moduli_buf[12], moduli_buf[13], moduli_buf[14]
 #ifndef UNIT16	/* and not UNIT8 */
-	,&moduli_buf[15][0]
-	,&moduli_buf[16][0], &moduli_buf[17][0], &moduli_buf[18][0], &moduli_buf[19][0], 
-	 &moduli_buf[20][0], &moduli_buf[21][0], &moduli_buf[22][0], &moduli_buf[23][0], 
-	 &moduli_buf[24][0], &moduli_buf[25][0], &moduli_buf[26][0], &moduli_buf[27][0], 
-	 &moduli_buf[28][0], &moduli_buf[29][0], &moduli_buf[30][0]
+	,moduli_buf[15] ,moduli_buf[16], moduli_buf[17], moduli_buf[18],
+	 moduli_buf[19], moduli_buf[20], moduli_buf[21], moduli_buf[22],
+	 moduli_buf[23], moduli_buf[24], moduli_buf[25], moduli_buf[26],
+	 moduli_buf[27], moduli_buf[28], moduli_buf[29], moduli_buf[30]
 #endif	/* UNIT16 and UNIT8 not defined */
 #endif	/* UNIT8 not defined */
 };
 
-/*	To optimize msubs, need following 2 unit arrays, each filled
-	with the most significant unit and next-to-most significant unit
-	of the preshifted images of the modulus. */
+/*
+ * To optimize msubs, need following 2 unit arrays, each filled
+ * with the most significant unit and next-to-most significant unit
+ * of the preshifted images of the modulus.
+ */
 static unit msu_moduli[UNITSIZE] = {0}; /* most signif. unit */
 static unit nmsu_moduli[UNITSIZE] = {0}; /* next-most signif. unit */
 
-/*	mpdbuf contains preshifted images of the multiplicand, mod n.
-	It is used only by mp_modmult.  It could be staticly declared
-	inside of mp_modmult, but we put it outside mp_modmult so that
-	it can be wiped clean by modmult_burn(), which is called at the
-	end of mp_modexp.  This is so that no sensitive data is left in
-	memory after the program exits.
-*/
+/*
+ * mpdbuf contains preshifted images of the multiplicand, mod n.
+ * It is used only by mp_modmult.  It could be staticly declared
+ * inside of mp_modmult, but we put it outside mp_modmult so that
+ * it can be wiped clean by modmult_burn(), which is called at the
+ * end of mp_modexp.  This is so that no sensitive data is left in
+ * memory after the program exits.
+ */
 static unit mpdbuf[UNITSIZE-1][MAX_UNIT_PRECISION] = {0};
 
 
 static void stage_mp_images(unitptr images[UNITSIZE],unitptr r)
-/*	Computes UNITSIZE images of r, each shifted left 1 more bit.
-	Used only by modmult function.
-*/
-{	short int i;
+/*
+ * Computes UNITSIZE images of r, each shifted left 1 more bit.
+ * Used only by modmult function.
+ */
+{
+	short int i;
+
 	images[0] = r;	/* no need to move the first image, just copy ptr */
-	for (i=1; i<UNITSIZE; i++)
-	{	mp_move(images[i],images[i-1]);
+	for (i=1; i<UNITSIZE; i++) {
+		mp_move(images[i],images[i-1]);
 		mp_shift_left(images[i]);
 		msub(images[i],moduli[0]);
 	}
-}	/* stage_mp_images */
+} /* stage_mp_images */
 
 
 int stage_merritt_modulus(unitptr n)
-/*	Computes UNITSIZE images of modulus, each shifted left 1 more bit.
-	Before calling mp_modmult, you must first stage the moduli images by
-	calling stage_modulus.  n is the pointer to the modulus.
-	Assumes that global_precision has already been adjusted to the
-	size of the modulus, plus SLOP_BITS.
-*/
-{	short int i;
+/*
+ * Computes UNITSIZE images of modulus, each shifted left 1 more bit.
+ * Before calling mp_modmult, you must first stage the moduli images by
+ * calling stage_modulus.  n is the pointer to the modulus.
+ * Assumes that global_precision has already been adjusted to the
+ * size of the modulus, plus SLOP_BITS.
+ */
+{
+	short int i;
 	unitptr msu;	/* ptr to most significant unit, for faster msubs */
+
 	moduli[0] = n;	/* no need to move the first image, just copy ptr */
 
 	/* used by optimized msubs macro... */
@@ -809,8 +847,8 @@ int stage_merritt_modulus(unitptr n)
 	msu_moduli[0] = *post_lowerunit(msu);
 	nmsu_moduli[0] = *msu;
 
-	for (i=1; i<UNITSIZE; i++)
-	{	mp_move(moduli[i],moduli[i-1]);
+	for (i=1; i<UNITSIZE; i++) {
+		mp_move(moduli[i],moduli[i-1]);
 		mp_shift_left(moduli[i]);
 
 		/* used by optimized msubs macro... */
@@ -818,8 +856,8 @@ int stage_merritt_modulus(unitptr n)
 		msu_moduli[i] = *post_lowerunit(msu);	/* for faster msubs */
 		nmsu_moduli[i] = *msu;
 	}
-	return(0);	/* normal return */
-}	/* stage_merritt_modulus */
+	return 0;	/* normal return */
+} /* stage_merritt_modulus */
 
 
 /* The following macros, sniffadd and msubs, are used by modmult... */
@@ -829,13 +867,14 @@ int stage_merritt_modulus(unitptr n)
 /* #define msubs0(i) msub(prod,moduli[i])
 */
 
-/*	To optimize msubs, compare the most significant units of the
-	product and the shifted modulus before deciding to call the full
-	mp_compare routine.  Better still, compare the upper two units
-	before deciding to call mp_compare.
-	Optimization of msubs relies heavily on standard C left-to-right
-	parsimonious evaluation of logical expressions.
-*/
+/*
+ * To optimize msubs, compare the most significant units of the
+ * product and the shifted modulus before deciding to call the full
+ * mp_compare routine.  Better still, compare the upper two units
+ * before deciding to call mp_compare.
+ * Optimization of msubs relies heavily on standard C left-to-right
+ * parsimonious evaluation of logical expressions.
+ */
 
 /* Partially-optimized msubs macro (msubs1) follows... */
 /* #define msubs1(i) if ( \
@@ -853,12 +892,13 @@ int stage_merritt_modulus(unitptr n)
 
 int merritt_modmult(register unitptr prod,
 	unitptr multiplicand,register unitptr multiplier)
-	/*	Performs combined multiply/modulo operation.
-		Computes:  prod = (multiplicand*multiplier) mod modulus
-		WARNING: All the arguments must be less than the modulus!
-		Assumes the modulus has been predefined by first calling
-		stage_modulus.
-	*/
+/*
+ * Performs combined multiply/modulo operation.
+ * Computes:  prod = (multiplicand*multiplier) mod modulus
+ * WARNING: All the arguments must be less than the modulus!
+ * Assumes the modulus has been predefined by first calling
+ * stage_modulus.
+ */
 {
 	/* p_m, msu_prod, and nmsu_prod are used by the optimized msubs macro...*/
 	register signedunit p_m;
@@ -867,17 +907,17 @@ int merritt_modmult(register unitptr pro
 	short mprec;		/* precision of multiplier, in units */
 	/*	Array mpd contains a list of pointers to preshifted images of
 		the multiplicand: */
-	static unitptr mpd[UNITSIZE] =
-	{	 0,              &mpdbuf[ 0][0], &mpdbuf[ 1][0], &mpdbuf[ 2][0],
-		 &mpdbuf[ 3][0], &mpdbuf[ 4][0], &mpdbuf[ 5][0], &mpdbuf[ 6][0]
+	static unitptr mpd[UNITSIZE] = {
+		 0,          mpdbuf[ 0], mpdbuf[ 1], mpdbuf[ 2],
+		 mpdbuf[ 3], mpdbuf[ 4], mpdbuf[ 5], mpdbuf[ 6]
 #ifndef UNIT8
-		,&mpdbuf[ 7][0], &mpdbuf[ 8][0], &mpdbuf[ 9][0], &mpdbuf[10][0],
-		 &mpdbuf[11][0], &mpdbuf[12][0], &mpdbuf[13][0], &mpdbuf[14][0]
+		,mpdbuf[ 7], mpdbuf[ 8], mpdbuf[ 9], mpdbuf[10],
+		 mpdbuf[11], mpdbuf[12], mpdbuf[13], mpdbuf[14]
 #ifndef UNIT16	/* and not UNIT8 */
-		,&mpdbuf[15][0], &mpdbuf[16][0], &mpdbuf[17][0], &mpdbuf[18][0],
-		 &mpdbuf[19][0], &mpdbuf[20][0], &mpdbuf[21][0], &mpdbuf[22][0],
-		 &mpdbuf[23][0], &mpdbuf[24][0], &mpdbuf[25][0], &mpdbuf[26][0],
-		 &mpdbuf[27][0], &mpdbuf[28][0], &mpdbuf[29][0], &mpdbuf[30][0]
+		,mpdbuf[15], mpdbuf[16], mpdbuf[17], mpdbuf[18],
+		 mpdbuf[19], mpdbuf[20], mpdbuf[21], mpdbuf[22],
+		 mpdbuf[23], mpdbuf[24], mpdbuf[25], mpdbuf[26],
+		 mpdbuf[27], mpdbuf[28], mpdbuf[29], mpdbuf[30]
 #endif	/* UNIT16 and UNIT8 not defined */
 #endif	/* UNIT8 not defined */
 	};
@@ -890,51 +930,54 @@ int merritt_modmult(register unitptr pro
 	nmsu_prod = msu_prod;
 	post_lowerunit(nmsu_prod); /* Get next-MSU of prod */
 
-	/*	To understand this algorithm, it would be helpful to first
-		study the conventional Russian peasant modmult algorithm.
-		This one does about the same thing as Russian peasant, but
-		is organized differently to save some steps.  It loops
-		through the multiplier a word (unit) at a time, instead of
-		a bit at a time.  It word-shifts the product instead of
-		bit-shifting it, so it should be faster.  It also does about
-		half as many subtracts as Russian peasant.
-	*/
+	/*
+	 * To understand this algorithm, it would be helpful to first
+	 * study the conventional Russian peasant modmult algorithm.
+	 * This one does about the same thing as Russian peasant, but
+	 * is organized differently to save some steps.  It loops
+	 * through the multiplier a word (unit) at a time, instead of
+	 * a bit at a time.  It word-shifts the product instead of
+	 * bit-shifting it, so it should be faster.  It also does about
+	 * half as many subtracts as Russian peasant.
+	 */
 
 	mp_init(prod,0);	/* Initialize product to 0. */
 
-	/*	The way mp_modmult is actually used in cryptographic
-		applications, there will NEVER be a zero multiplier or
-		multiplicand.  So there is no need to optimize for that
-		condition.
-	*/
+	/*
+	 * The way mp_modmult is actually used in cryptographic
+	 * applications, there will NEVER be a zero multiplier or
+	 * multiplicand.  So there is no need to optimize for that
+	 * condition.
+	 */
 /*	if (testeq(multiplicand,0))
-		return(0); */	/* zero multiplicand means zero product */
+		return 0; */	/* zero multiplicand means zero product */
 	/* Normalize and compute number of units in multiplier first: */
 	normalize(multiplier,mprec);
 	if (mprec==0)	/* if precision of multiplier is 0 */
-		return(0);	/* zero multiplier means zero product */
+		return 0;	/* zero multiplier means zero product */
 	make_msbptr(multiplier,mprec); /* start at MSU of multiplier */
 
-	while (mprec--)	/* Loop for the number of units in the multiplier */
-	{
-		/*	Shift the product one whole unit to the left.
-			This is part of the multiply phase of modmult.
-		*/
+	while (mprec--) { /* Loop for the number of units in the multiplier */
+		/*
+		 * Shift the product one whole unit to the left.
+		 * This is part of the multiply phase of modmult.
+		 */
 
 		mp_lshift_unit(prod);
 
-		/*      The product may have grown by as many as UNITSIZE
-			bits.  That's why we have global_precision set to the
-			size of the modulus plus UNITSIZE slop bits.
-			Now reduce the product back down by conditionally
-			subtracting the preshifted images of the modulus.
-			This is part of the modulo reduction phase of modmult. 
-			The following loop is unrolled for speed, using macros...
+		/* 
+		 * The product may have grown by as many as UNITSIZE
+		 * bits.  That's why we have global_precision set to the
+		 * size of the modulus plus UNITSIZE slop bits.
+		 * Now reduce the product back down by conditionally
+		 * subtracting the preshifted images of the modulus.
+		 * This is part of the modulo reduction phase of modmult. 
+		 * The following loop is unrolled for speed, using macros...
 
 		for (i=UNITSIZE-1; i>=LOG_UNITSIZE; i--)
 			if (mp_compare(prod,moduli[i]) >= 0)
 				mp_sub(prod,moduli[i]);
-		*/
+		 */
 
 #ifndef UNIT8
 #ifndef UNIT16  /* and not UNIT8 */
@@ -974,17 +1017,18 @@ int merritt_modmult(register unitptr pro
 #endif
 #endif
 
-		/*	Sniff each bit in the current unit of the multiplier, 
-			and conditionally add the corresponding preshifted 
-			image of the multiplicand to the product.
-			This is also part of the multiply phase of modmult.
-
-			The following loop is unrolled for speed, using macros...
+		/*
+		 * Sniff each bit in the current unit of the multiplier, 
+		 * and conditionally add the corresponding preshifted 
+		 * image of the multiplicand to the product.
+		 * This is also part of the multiply phase of modmult.
+		 *
+		 * The following loop is unrolled for speed, using macros...
 
 		for (i=UNITSIZE-1; i>=0; i--)
 		   if (*multiplier & power_of_2(i))
 				mp_add(prod,mpd[i]);
-		*/
+		 */
 #ifndef UNIT8
 #ifndef UNIT16	/* and not UNIT8 */
 		sniffadd(31);
@@ -1022,18 +1066,19 @@ int merritt_modmult(register unitptr pro
 		sniffadd(1);
 		sniffadd(0);
 
-		/*	The product may have grown by as many as LOG_UNITSIZE+1
-			bits.
-			Now reduce the product back down by conditionally 
-			subtracting LOG_UNITSIZE+1 preshifted images of the 
-			modulus.  This is the modulo reduction phase of modmult.
-
-			The following loop is unrolled for speed, using macros...
+		/*
+		 * The product may have grown by as many as LOG_UNITSIZE+1
+		 * bits.
+		 * Now reduce the product back down by conditionally 
+		 * subtracting LOG_UNITSIZE+1 preshifted images of the 
+		 * modulus.  This is the modulo reduction phase of modmult.
+		 *
+		 * The following loop is unrolled for speed, using macros...
 
 		for (i=LOG_UNITSIZE; i>=0; i--)
 			if (mp_compare(prod,moduli[i]) >= 0) 
 				mp_sub(prod,moduli[i]); 
-		*/
+		 */
 
 #ifndef UNIT8
 #ifndef UNIT16
@@ -1049,25 +1094,27 @@ int merritt_modmult(register unitptr pro
 		/* Bump pointer to next lower unit of multiplier: */
 		post_lowerunit(multiplier);
 
-	}	/* Loop for the number of units in the multiplier */
+	} /* Loop for the number of units in the multiplier */
 
-	return(0);	/* normal return */
+	return 0;	/* normal return */
 
-}	/* merritt_modmult */
+} /* merritt_modmult */
 
 
 #undef msubs
 #undef sniffadd
 
 
-/*	Merritt's mp_modmult function leaves some internal tables in memory,
-	so we have to call modmult_burn() at the end of mp_modexp.
-	This is so that no cryptographically sensitive data is left in memory
-	after the program exits.
-*/
+/*
+ * Merritt's mp_modmult function leaves some internal tables in memory,
+ * so we have to call modmult_burn() at the end of mp_modexp.
+ * This is so that no cryptographically sensitive data is left in memory
+ * after the program exits.
+ */
 void merritt_burn(void)
-/*	Alias for modmult_burn, merritt_burn() is called only by mp_modexp. */
-{	unitfill0(&(mpdbuf[0][0]),(UNITSIZE-1)*MAX_UNIT_PRECISION);
+/* Alias for modmult_burn, merritt_burn() is called only by mp_modexp. */
+{
+	unitfill0(&(mpdbuf[0][0]),(UNITSIZE-1)*MAX_UNIT_PRECISION);
 	unitfill0(&(moduli_buf[0][0]),(UNITSIZE-1)*MAX_UNIT_PRECISION);
 	unitfill0(msu_moduli,UNITSIZE);
 	unitfill0(nmsu_moduli,UNITSIZE);
@@ -1080,38 +1127,42 @@ void merritt_burn(void)
 
 #ifdef UPTON_OR_SMITH	/* Used by Upton's and Smith's modmult algorithms */
 
-/*	Jimmy Upton's multiprecision modmult algorithm in C.
-	Performs a multiply combined with a modulo operation.
-
-	The following support functions and data structures
-	are used only by Upton's modmult algorithm...
-*/
+/*
+ * Jimmy Upton's multiprecision modmult algorithm in C.
+ * Performs a multiply combined with a modulo operation.
+ *
+ * The following support functions and data structures
+ * are used only by Upton's modmult algorithm...
+ */
 
 short munit_prec;	/* global_precision expressed in MULTUNITs */
 
-/*	Note that since the SPARC CPU has no hardware integer multiply
-	instruction, there is not that much advantage in having an
-	assembly version of mp_smul on that machine.  It might be faster
-	to use Merritt's modmult instead of Upton's modmult on the SPARC.
-*/
+/*
+ * Note that since the SPARC CPU has no hardware integer multiply
+ * instruction, there is not that much advantage in having an
+ * assembly version of mp_smul on that machine.  It might be faster
+ * to use Merritt's modmult instead of Upton's modmult on the SPARC.
+ */
 
 /*
-	Multiply the single-word multiplier times the multiprecision integer
-	in multiplicand, accumulating result in prod.  The resulting
-	multiprecision prod will be 1 word longer than the multiplicand.
-	multiplicand is munit_prec words long.  We add into prod, so caller
-	should zero it out first.  For best results, this time-critical
-	function should be implemented in assembly.
-	NOTE:  Unlike other functions in the multiprecision arithmetic
-	library, both multiplicand and prod are pointing at the LSB,
-	regardless of byte order of the machine.  On an 80x86, this makes
-	no difference.  But if this assembly function is implemented
-	on a 680x0, it becomes important.
-*/
-/*	Note that this has been modified from the previous version to allow
-	better support for Smith's modmult:
-		The final carry bit is added to the existing product
-		array, rather than simply stored.
+ * Multiply the single-word multiplier times the multiprecision integer
+ * in multiplicand, accumulating result in prod.  The resulting
+ * multiprecision prod will be 1 word longer than the multiplicand.
+ * multiplicand is munit_prec words long.  We add into prod, so caller
+ * should zero it out first.  For best results, this time-critical
+ * function should be implemented in assembly.
+ * NOTE:  Unlike other functions in the multiprecision arithmetic
+ * library, both multiplicand and prod are pointing at the LSB,
+ * regardless of byte order of the machine.  On an 80x86, this makes
+ * no difference.  But if this assembly function is implemented
+ * on a 680x0, it becomes important.
+ */
+
+/*
+ * Note that this has been modified from the previous version to allow
+ * better support for Smith's modmult:
+ * 	The final carry bit is added to the existing product
+ * 	array, rather than simply stored.
 */
 
 #ifndef mp_smul
@@ -1121,22 +1172,23 @@ void mp_smul (MULTUNIT *prod, MULTUNIT *
 	unsigned long p, carry;
 
 	carry = 0;
-	for (i=0; i<munit_prec; ++i)
-	{	p = (unsigned long)multiplier * *post_higherunit(multiplicand);
+	for (i=0; i<munit_prec; ++i) {
+		p = (unsigned long)multiplier * *post_higherunit(multiplicand);
 		p += *prod + carry;
 		*post_higherunit(prod) = (MULTUNIT)p;
 		carry = p >> MULTUNITSIZE;
 	}
 	/* Add carry to the next higher word of product / dividend */
 	*prod += (MULTUNIT)carry;
-}	/* mp_smul */
+} /* mp_smul */
 #endif
 
-/*	mp_dmul is a double-precision multiply multiplicand times multiplier,
-	result into prod.  prod must be pointing at a "double precision"
-	buffer.  E.g. If global_precision is 10 words, prod must be
-	pointing at a 20-word buffer.
-*/
+/*
+ * mp_dmul is a double-precision multiply multiplicand times multiplier,
+ * result into prod.  prod must be pointing at a "double precision"
+ * buffer.  E.g. If global_precision is 10 words, prod must be
+ * pointing at a 20-word buffer.
+ */
 #ifndef mp_dmul
 void mp_dmul (unitptr prod, unitptr multiplicand, unitptr multiplier)
 {
@@ -1158,7 +1210,7 @@ void mp_dmul (unitptr prod, unitptr mult
 	for (i=0; i<munit_prec; ++i)
 		mp_smul (post_higherunit(prodp),
 			p_multiplicand, *post_higherunit(p_multiplier));
-}	/* mp_dmul */
+} /* mp_dmul */
 #endif  /* mp_dmul */
 
 static unit ALIGN modulus[MAX_UNIT_PRECISION];
@@ -1167,12 +1219,13 @@ static short nbits;			/* number of modul
 
 #ifdef UPTON
 
-/*	These scratchpad arrays are used only by upton_modmult (mp_modmult).
-	Some of them could be staticly declared inside of mp_modmult, but we
-	put them outside mp_modmult so that they can be wiped clean by
-	modmult_burn(), which is called at the end of mp_modexp.  This is
-	so that no sensitive data is left in memory after the program exits.
-*/
+/*
+ * These scratchpad arrays are used only by upton_modmult (mp_modmult).
+ * Some of them could be staticly declared inside of mp_modmult, but we
+ * put them outside mp_modmult so that they can be wiped clean by
+ * modmult_burn(), which is called at the end of mp_modexp.  This is
+ * so that no sensitive data is left in memory after the program exits.
+ */
 
 static unit ALIGN reciprocal[MAX_UNIT_PRECISION];
 static unit ALIGN dhi[MAX_UNIT_PRECISION];
@@ -1183,13 +1236,14 @@ static unit ALIGN f_data[MAX_UNIT_PRECIS
 static short nbitsDivUNITSIZE;
 static short nbitsModUNITSIZE;
 
-/*	stage_upton_modulus() is aliased to stage_modulus().
-	Prepare for an Upton modmult.  Calculate the reciprocal of modulus,
-	and save both.  Note that reciprocal will have one more bit than
-	modulus.
-	Assumes that global_precision has already been adjusted to the
-	size of the modulus, plus SLOP_BITS.
-*/
+/*
+ * stage_upton_modulus() is aliased to stage_modulus().
+ * Prepare for an Upton modmult.  Calculate the reciprocal of modulus,
+ * and save both.  Note that reciprocal will have one more bit than
+ * modulus.
+ * Assumes that global_precision has already been adjusted to the
+ * size of the modulus, plus SLOP_BITS.
+ */
 int stage_upton_modulus(unitptr n)
 {
 	mp_move(modulus, n);
@@ -1197,17 +1251,18 @@ int stage_upton_modulus(unitptr n)
 	nbits = countbits(modulus);
 	nbitsDivUNITSIZE = nbits / UNITSIZE;
 	nbitsModUNITSIZE = nbits % UNITSIZE;
-	return(0);	/* normal return */
-}	/* stage_upton_modulus */
+	return 0;	/* normal return */
+} /* stage_upton_modulus */
 
 
-/*	Upton's algorithm performs a multiply combined with a modulo operation.
-	Computes:  prod = (multiplicand*multiplier) mod modulus
-	WARNING: All the arguments must be less than the modulus!
-	References global unitptr modulus and reciprocal.
-	The reciprocal of modulus is 1 bit longer than the modulus.
-	upton_modmult() is aliased to mp_modmult().
-*/
+/*
+ * Upton's algorithm performs a multiply combined with a modulo operation.
+ * Computes:  prod = (multiplicand*multiplier) mod modulus
+ * WARNING: All the arguments must be less than the modulus!
+ * References global unitptr modulus and reciprocal.
+ * The reciprocal of modulus is 1 bit longer than the modulus.
+ * upton_modmult() is aliased to mp_modmult().
+ */
 int upton_modmult (unitptr prod, unitptr multiplicand, unitptr multiplier)
 {
 	unitptr d = d_data;
@@ -1252,8 +1307,8 @@ int upton_modmult (unitptr prod, unitptr
 		mp_sub (d,modulus);
 
 	mp_move(prod,d);
-	return(0);  /* normal return */
-}	/* upton_modmult */
+	return 0;  /* normal return */
+} /* upton_modmult */
 
 
 /*	Upton's mp_modmult function leaves some internal arrays in memory,
@@ -1271,7 +1326,7 @@ void upton_burn(void)
 	unitfill0(e_data,MAX_UNIT_PRECISION*2);
 	unitfill0(f_data,MAX_UNIT_PRECISION*2);
 	nbits = nbitsDivUNITSIZE = nbitsModUNITSIZE = 0;
-}	/* upton_burn */
+} /* upton_burn */
 
 /******* end of Upton's MODMULT stuff. *******/
 /*=========================================================================*/
@@ -1279,27 +1334,29 @@ void upton_burn(void)
 
 #ifdef SMITH	/* using Thad Smith's modmult algorithm */
 
-/*	Thad Smith's implementation of multiprecision modmult algorithm in C.
-	Performs a multiply combined with a modulo operation.
-	The multiplication is done with mp_dmul, the same as for Upton's
-	modmult.  The modulus reduction is done by long division, in
-	which a trial quotient "digit" is determined, then the product of
-	that digit and the divisor are subtracted from the dividend.
-
-	In this case, the digit is MULTUNIT in size and the subtraction
-	is done by adding the product to the one's complement of the
-	dividend, which allows use of the existing mp_smul routine.
-
-	The following support functions and data structures
-	are used only by Smith's modmult algorithm...
-*/
+/*
+ * Thad Smith's implementation of multiprecision modmult algorithm in C.
+ * Performs a multiply combined with a modulo operation.
+ * The multiplication is done with mp_dmul, the same as for Upton's
+ * modmult.  The modulus reduction is done by long division, in
+ * which a trial quotient "digit" is determined, then the product of
+ * that digit and the divisor are subtracted from the dividend.
+ *
+ * In this case, the digit is MULTUNIT in size and the subtraction
+ * is done by adding the product to the one's complement of the
+ * dividend, which allows use of the existing mp_smul routine.
+ *
+ * The following support functions and data structures
+ * are used only by Smith's modmult algorithm...
+ */
 
-/*	These scratchpad arrays are used only by smith_modmult (mp_modmult).
-	Some of them could be statically declared inside of mp_modmult, but we
-	put them outside mp_modmult so that they can be wiped clean by
-	modmult_burn(), which is called at the end of mp_modexp.  This is
-	so that no sensitive data is left in memory after the program exits.
-*/
+/*
+ * These scratchpad arrays are used only by smith_modmult (mp_modmult).
+ * Some of them could be statically declared inside of mp_modmult, but we
+ * put them outside mp_modmult so that they can be wiped clean by
+ * modmult_burn(), which is called at the end of mp_modexp.  This is
+ * so that no sensitive data is left in memory after the program exits.
+ */
 
 static unit ALIGN ds_data[MAX_UNIT_PRECISION*2+2];
 
@@ -1316,12 +1373,13 @@ static MULTUNIT recipl; 		/* LSunit of s
 static short modlenMULTUNITS;		/* length of modulus in MULTUNITs */
 static MULTUNIT mutemp; 		/* temporary */
 
-/*	The routines mp_smul and mp_dmul are the same as for UPTON and
-	should be coded in assembly.  Note, however, that the previous
-	Upton's mp_smul version has been modified to compatible with
-	Smith's modmult.  The new version also still works for Upton's
-	modmult.
-*/
+/*
+ * The routines mp_smul and mp_dmul are the same as for UPTON and
+ * should be coded in assembly.  Note, however, that the previous
+ * Upton's mp_smul version has been modified to compatible with
+ * Smith's modmult.  The new version also still works for Upton's
+ * modmult.
+ */
 
 #ifndef mp_set_recip
 #define mp_set_recip(rh,rl,m)     /* null */
@@ -1335,71 +1393,75 @@ MULTUNIT mp_quo_digit (MULTUNIT *dividen
 #define mp_musubb mp_subb	/* use existing routine */
 #else  /* ! MULTUNIT_SIZE_SAME */
 
-/*	This performs the same function as mp_subb, but with MULTUNITs.
-	Note: Processors without alignment requirements may be able to use
-	mp_subb, even though MULTUNITs are smaller than units.  In that case,
-	use mp_subb, since it would be faster if coded in assembly.  Note that
-	this implementation won't work for MULTUNITs which are long -- use
-	mp_subb in that case.
-*/
+/*
+ * This performs the same function as mp_subb, but with MULTUNITs.
+ * Note: Processors without alignment requirements may be able to use
+ * mp_subb, even though MULTUNITs are smaller than units.  In that case,
+ * use mp_subb, since it would be faster if coded in assembly.  Note that
+ * this implementation won't work for MULTUNITs which are long -- use
+ * mp_subb in that case.
+ */
 #ifndef mp_musubb
 boolean mp_musubb
 	(register MULTUNIT* r1,register MULTUNIT* r2,register boolean borrow)
-	/* multiprecision subtract of MULTUNITs with borrow, r2 from r1,
-	** result in r1 */
-	/* borrow is incoming borrow flag-- value should be 0 or 1 */
-{	register ulint x;	/* won't work if sizeof(MULTUNIT)==
+/*
+ * multiprecision subtract of MULTUNITs with borrow, r2 from r1, result in r1
+ * borrow is incoming borrow flag-- value should be 0 or 1
+ */
+{
+	register ulint x;	/* won't work if sizeof(MULTUNIT)==
 						 sizeof(long)	    */
 	short precision;	/* number of MULTUNITs to subtract */
 	precision = global_precision * MULTUNITs_perunit;
 	make_lsbptr(r1,precision);
 	make_lsbptr(r2,precision);
-	while (precision--)
-	{	x = (ulint) *r1 - (ulint) *post_higherunit(r2) - (ulint) borrow;
+	while (precision--) {
+		x = (ulint) *r1 - (ulint) *post_higherunit(r2) - (ulint) borrow;
 		*post_higherunit(r1) = x;
 		borrow = (((1L << MULTUNITSIZE) & x) != 0L);
 	}
 	return (borrow);
-}	/* mp_musubb */
+} /* mp_musubb */
 #endif	/* mp_musubb */
 #endif	/* MULTUNIT_SIZE_SAME */
 
-/*	The function mp_quo_digit is the heart of Smith's modulo reduction,
-	which uses a form of long division.  It computes a trial quotient
-	"digit" (MULTUNIT-sized digit) by multiplying the three most
-	significant MULTUNITs of the dividend by the two most significant
-	MULTUNITs of the reciprocal of the modulus.  Note that this function
-	requires that MULTUNITSIZE * 2 <= sizeof(unsigned long).
-
-	An important part of this technique is that the quotient never be
-	too small, although it may occasionally be too large.  This was
-	done to eliminate the need to check and correct for a remainder
-	exceeding the divisor.	It is easier to check for a negative
-	remainder.  The following technique rarely needs correction for
-	MULTUNITs of at least 16 bits.
-
-	The following routine has two implementations:
-
-	ASM_PROTOTYPE defined: written to be an executable prototype for
-	an efficient assembly language implementation.  Note that several
-	of the following masks and shifts can be done by directly
-	manipulating single precision registers on some architectures.
-
-	ASM_PROTOTYPE undefined: a slightly more efficient implementation
-	in C.  Although this version returns a result larger than the
-	optimum (which is corrected in smith_modmult) more often than the
-	prototype version, the faster execution in C more than makes up
-	for the difference.
- 
-	Parameter: dividend - points to the most significant MULTUNIT
-		of the dividend.  Note that dividend actually contains the 
-		one's complement of the actual dividend value (see comments for 
-		smith_modmult).
-	
-	 Return: the trial quotient digit resulting from dividing the first
-		three MULTUNITs at dividend by the upper two MULTUNITs of the 
-		modulus.
-*/
+/*
+ * The function mp_quo_digit is the heart of Smith's modulo reduction,
+ * which uses a form of long division.  It computes a trial quotient
+ * "digit" (MULTUNIT-sized digit) by multiplying the three most
+ * significant MULTUNITs of the dividend by the two most significant
+ * MULTUNITs of the reciprocal of the modulus.  Note that this function
+ * requires that MULTUNITSIZE * 2 <= sizeof(unsigned long).
+ *
+ * An important part of this technique is that the quotient never be
+ * too small, although it may occasionally be too large.  This was
+ * done to eliminate the need to check and correct for a remainder
+ * exceeding the divisor.	It is easier to check for a negative
+ * remainder.  The following technique rarely needs correction for
+ * MULTUNITs of at least 16 bits.
+ *
+ * The following routine has two implementations:
+ *
+ * ASM_PROTOTYPE defined: written to be an executable prototype for
+ * an efficient assembly language implementation.  Note that several
+ * of the following masks and shifts can be done by directly
+ * manipulating single precision registers on some architectures.
+ *
+ * ASM_PROTOTYPE undefined: a slightly more efficient implementation
+ * in C.  Although this version returns a result larger than the
+ * optimum (which is corrected in smith_modmult) more often than the
+ * prototype version, the faster execution in C more than makes up
+ * for the difference.
+ * 
+ * Parameter: dividend - points to the most significant MULTUNIT
+ * 	of the dividend.  Note that dividend actually contains the 
+ * 	one's complement of the actual dividend value (see comments for 
+ * 	smith_modmult).
+ * 
+ *  Return: the trial quotient digit resulting from dividing the first
+ * 	three MULTUNITs at dividend by the upper two MULTUNITs of the 
+ * 	modulus.
+ */
 
 /* #define ASM_PROTOTYPE */ /* undefined: use C-optimized version */
 
@@ -1408,10 +1470,11 @@ MULTUNIT mp_quo_digit (MULTUNIT *dividen
 	unsigned long q, q0, q1, q2;
 	unsigned short lsb_factor;
 
-/*	Compute the least significant product group.
-	The last terms of q1 and q2 perform upward rounding, which is
-	needed to guarantee that the result not be too small.
-*/
+/*
+ * Compute the least significant product group.
+ *The last terms of q1 and q2 perform upward rounding, which is
+ *needed to guarantee that the result not be too small.
+ */
 	q1 = (dividend[tohigher(-2)] ^ MULTUNIT_mask) * (unsigned long)reciph
 	     + reciph;
 	q2 = (dividend[tohigher(-1)] ^ MULTUNIT_mask) * (unsigned long)recipl
@@ -1459,19 +1522,20 @@ MULTUNIT mp_quo_digit (MULTUNIT *dividen
 }
 #endif	/* mp_quo_digit */
 
-/*	stage_smith_modulus() - Prepare for a Smith modmult.
-
-	Calculate the reciprocal of modulus with a precision of two MULTUNITs.
-	Assumes that global_precision has already been adjusted to the
-	size of the modulus, plus SLOP_BITS.
-
-	Note: This routine was designed to work with large values and
-	doesn't have the necessary testing or handling to work with a
-	modulus having less than three significant units.  For such cases,
-	the separate multiply and modulus routines can be used.
-
-	stage_smith_modulus() is aliased to stage_modulus().
-*/
+/*
+ * stage_smith_modulus() - Prepare for a Smith modmult.
+ * 
+ * Calculate the reciprocal of modulus with a precision of two MULTUNITs.
+ * Assumes that global_precision has already been adjusted to the
+ * size of the modulus, plus SLOP_BITS.
+ * 
+ * Note: This routine was designed to work with large values and
+ * doesn't have the necessary testing or handling to work with a
+ * modulus having less than three significant units.  For such cases,
+ * the separate multiply and modulus routines can be used.
+ * 
+ * stage_smith_modulus() is aliased to stage_modulus().
+ */
 
 int stage_smith_modulus(unitptr n_modulus)
 {
@@ -1527,8 +1591,8 @@ int stage_smith_modulus(unitptr n_modulu
 	reciph = *mpm;
 	mp_set_recip (reciph, recipl, mshift);
 	set_precision (original_precision);
-	return(0);	/* normal return */
-}	/* stage_smith_modulus */
+	return 0;	/* normal return */
+} /* stage_smith_modulus */
 
 /*	Smith's algorithm performs a multiply combined with a modulo operation.
 	Computes:  prod = (multiplicand*multiplier) mod modulus
@@ -1564,8 +1628,8 @@ smith_modmult(unitptr prod, unitptr mult
 	have to process leading zeros in our divide loop.
 */
 	dmi = significance(d) * MULTUNITs_perunit;
-	if (dmi >= modlenMULTUNITS)
-	{	/* Make dividend negative.  This allows the use of mp_smul to
+	if (dmi >= modlenMULTUNITS) {
+		/* Make dividend negative.  This allows the use of mp_smul to
 		 * "subtract" the product of the modulus and the trial divisor
 		 * by actually adding to a negative dividend.
 		 * The one's complement of the dividend is used, since it causes
@@ -1585,26 +1649,27 @@ smith_modmult(unitptr prod, unitptr mult
 		nqd = dmi + 1 - modlenMULTUNITS;
 		dmpl = dmph - tohigher(modlenMULTUNITS);
 
-/*	Divide loop.
-	Each iteration computes the next quotient MULTUNIT digit, then
-	multiplies the divisor (modulus) by the quotient digit and adds
-	it to the one's complement of the dividend (equivalent to
-	subtracting).  If the product was greater than the remaining dividend,
-	we get a non-negative result, in which case we subtract off the
-	modulus to get the proper negative remainder.
-*/
-		for (; nqd; nqd--)
-		{	MULTUNIT    q;	    /* quotient trial digit */
+/*	
+ * Divide loop.
+ * Each iteration computes the next quotient MULTUNIT digit, then
+ * multiplies the divisor (modulus) by the quotient digit and adds
+ * it to the one's complement of the dividend (equivalent to
+ * subtracting).  If the product was greater than the remaining dividend,
+ * we get a non-negative result, in which case we subtract off the
+ * modulus to get the proper negative remainder.
+ */
+		for (; nqd; nqd--) {
+			MULTUNIT    q;	    /* quotient trial digit */
 
 			q = mp_quo_digit(dmph);
-			if (q > 0)
-			{	mp_smul(dmpl, modmpl, q);
+			if (q > 0) {
+				mp_smul(dmpl, modmpl, q);
 
 				/* Perform correction if q too large.
 				*  This rarely occurs.
 				*/
-				if (!(*dmph & MULTUNIT_msb))
-				{	dmp = dmpl;
+				if (!(*dmph & MULTUNIT_msb)) {
+					dmp = dmpl;
 					unmake_lsbptr(dmp, orig_precision * 
 						   MULTUNITs_perunit);
 					if (mp_musubb(dmp,
@@ -1620,14 +1685,14 @@ smith_modmult(unitptr prod, unitptr mult
 		set_precision(orig_precision);
 		mp_neg(d);
 		mp_dec(d);
-	} else
-	{	/* Product was less than modulus.  Return it. */
+	} else {
+		/* Product was less than modulus.  Return it. */
 		rescale(d, orig_precision * 2 + 1, orig_precision);
 		set_precision(orig_precision);
 	}
 	mp_move(prod, d);
 	return (0);		/* normal return */
-}				/* smith_modmult */
+} /* smith_modmult */
 
 
 /*	Smith's mp_modmult function leaves some internal arrays in memory,
@@ -1647,7 +1712,7 @@ void smith_burn(void)
 	reciph = recipl = 0;
 	modlenMULTUNITS = mutemp = 0;
 	mp_set_recip (0,0,0);
-}	/* smith_burn */
+} /* smith_burn */
 
 /*	End of Thad Smith's implementation of modmult. */
 
@@ -1656,11 +1721,12 @@ void smith_burn(void)
 
 int countbits(unitptr r)
 	/* Returns number of significant bits in r */
-{	int bits;
+{
+	int bits;
 	short prec;
 	register unit bitmask;
 	init_bitsniffer(r,bitmask,prec,bits);
-	return(bits);
+	return bits;
 } /* countbits */
 
 
@@ -1671,11 +1737,13 @@ char *copyright_notice(void)
 
 int mp_modexp(register unitptr expout,register unitptr expin,
 	register unitptr exponent,register unitptr modulus)
-{	/*	Russian peasant combined exponentiation/modulo algorithm.
-		Calls modmult instead of mult.
-		Computes:  expout = (expin**exponent) mod modulus
-		WARNING: All the arguments must be less than the modulus!
-	*/
+/*
+ * Russian peasant combined exponentiation/modulo algorithm.
+ * Calls modmult instead of mult.
+ * Computes:  expout = (expin**exponent) mod modulus
+ * WARNING: All the arguments must be less than the modulus!
+ */
+{
 	int bits;
 	short oldprecision;
 	register unit bitmask;
@@ -1687,21 +1755,21 @@ int mp_modexp(register unitptr expout,re
 	tally_modsquares = 0;	/* clear "number of modsquares" counter */
 #endif	/* COUNTMULTS */
 	mp_init(expout,1);
-	if (testeq(exponent,0))
-	{	if (testeq(expin,0))
-			return(-1); /* 0 to the 0th power means return error */
-		return(0);	/* otherwise, zero exponent means expout is 1 */
+	if (testeq(exponent,0)) {
+		if (testeq(expin,0))
+			return -1; /* 0 to the 0th power means return error */
+		return 0;	/* otherwise, zero exponent means expout is 1 */
 	}
 	if (testeq(modulus,0))
-		return(-2);		/* zero modulus means error */
+		return -2;		/* zero modulus means error */
 #if SLOP_BITS > 0	/* if there's room for sign bits */
 	if (mp_tstminus(modulus))
-		return(-2);		/* negative modulus means error */
+		return -2;		/* negative modulus means error */
 #endif	/* SLOP_BITS > 0 */
 	if (mp_compare(expin,modulus) >= 0)
-		return(-3); /* if expin >= modulus, return error */
+		return -3; /* if expin >= modulus, return error */
 	if (mp_compare(exponent,modulus) >= 0)
-		return(-4); /* if exponent >= modulus, return error */
+		return -4; /* if exponent >= modulus, return error */
 
 	oldprecision = global_precision;	/* save global_precision */
 	/* set smallest optimum precision for this modulus */
@@ -1712,9 +1780,9 @@ int mp_modexp(register unitptr expout,re
 	rescale(exponent,oldprecision,global_precision);
 	rescale(expout,oldprecision,global_precision);
 
-	if (stage_modulus(modulus))
-	{	set_precision(oldprecision);	/* restore original precision */
-		return(-5);		/* unstageable modulus (STEWART algorithm) */
+	if (stage_modulus(modulus)) {
+		set_precision(oldprecision);	/* restore original precision */
+		return -5;		/* unstageable modulus (STEWART algorithm) */
 	}
 
 	/* normalize and compute number of bits in exponent first */
@@ -1725,29 +1793,28 @@ int mp_modexp(register unitptr expout,re
 	mp_move(expout,expin);		/*  expout = (1*1)*expin; */
 	bump_bitsniffer(exponent,bitmask);
 
-	while (bits--)
-	{
+	while (bits--) {
 		poll_for_break(); /* polls keyboard, allows ctrl-C to abort program */
 #ifdef COUNTMULTS
 		tally_modsquares++;	/* bump "number of modsquares" counter */
 #endif	/* COUNTMULTS */
 		mp_modsquare(product,expout);
-		if (sniff_bit(exponent,bitmask))
-		{	mp_modmult(expout,product,expin);
+		if (sniff_bit(exponent,bitmask)) {
+			mp_modmult(expout,product,expin);
 #ifdef COUNTMULTS
 			tally_modmults++;	/* bump "number of modmults" counter */
 #endif	/* COUNTMULTS */
-		} else
-		{
+		} else {
 			mp_move(expout,product);
 		}
 		bump_bitsniffer(exponent,bitmask);
-	}	/* while bits-- */
+	} /* while bits-- */
 	mp_burn(product);	/* burn the evidence on the stack */
 	modmult_burn(); /* ask mp_modmult to also burn its own evidence */
 
 #ifdef COUNTMULTS	/* diagnostic analysis */
-	{	long atomic_mults;
+	{
+		long atomic_mults;
 		unsigned int unitcount,totalmults;
 		unitcount = bits2units(countbits(modulus));
 		/* calculation assumes modsquare takes as long as a modmult: */
@@ -1767,62 +1834,58 @@ int mp_modexp(register unitptr expout,re
 	   will be forced to include it in the executable object image... */
 	copyright_notice();	/* has no real effect at run time */
 
-	return(0);		/* normal return */
-}	/* mp_modexp */
+	return 0;		/* normal return */
+} /* mp_modexp */
 
 int mp_modexp_crt(unitptr expout, unitptr expin,
 	unitptr p, unitptr q, unitptr ep, unitptr eq, unitptr u)
-	/*	This is a faster modexp for moduli with a known
-		factorisation into two relatively prime factors p and q,
-		and an input relatively prime to the modulus,
-		the Chinese Remainder Theorem to do the computation
-		mod p and mod q, and then combine the results.  This
-		relies on a number of precomputed values, but does not
-		actually require the modulus n or the exponent e.
-
-		expout = expin ^ e mod (p*q).
-		We form this by evaluating
-		p2 = (expin ^ e) mod p and
-		q2 = (expin ^ e) mod q
-		and then combining the two by the CRT.
-
-		Two optimisations of this are possible.  First, we can
-		reduce expin modulo p and q before starting.
-
-		Second, since we know the factorisation of p and q
-		(trivially derived from the factorisation of n = p*q),
-		and expin is relatively prime to both p and q,
-		we can use Euler's theorem, expin^phi(m) = 1 (mod m),
-		to throw away multiples of phi(p) or phi(q) in e.
-		Letting ep = e mod phi(p) and
-		        eq = e mod phi(q)
-		then combining these two speedups, we only need to evaluate
-		p2 = ((expin mod p) ^ ep) mod p and
-		q2 = ((expin mod q) ^ eq) mod q.
-
-		Now we need to apply the CRT.  Starting with
-		expout = p2 (mod p) and
-		expout = q2 (mod q)
-		we can say that expout = p2 + p * k, and if we assume that
-		0 <= p2 < p, then 0 <= expout < p*q for some 0 <= k < q.
-		Since we want expout = q2 (mod q), then
-		p*k = q2-p2 (mod q).  Since p and q are relatively
-		prime, p has a multiplicative inverse u mod q.  In other
-		words,
-		       u = 1/p (mod q).
-		Multiplying by u on both sides gives k = u*(q2-p2) (mod q).
-		Since we want 0 <= k < q, we can thus find k as
-		k = (u * (q2-p2)) mod q.
-
-		Once we have k, evaluating p2 + p * k is easy, and
-		that gives us the result.
-
-		In the detailed implementation, there is a temporary, temp,
-		used to hold intermediate results, p2 is held in expout,
-		and q2 is used as a temporary in the final step when it is
-		no longer needed.  With that, you should be able to
-		understand the code below.
-	*/
+/*
+ * This is a faster modexp for moduli with a known factorisation into two
+ * relatively prime factors p and q, and an input relatively prime to the
+ * modulus, the Chinese Remainder Theorem to do the computation mod p and
+ * mod q, and then combine the results.  This relies on a number of
+ * precomputed values, but does not actually require the modulus n or the
+ * exponent e.
+ * 
+ * expout = expin ^ e mod (p*q).
+ * We form this by evaluating
+ * p2 = (expin ^ e) mod p and
+ * q2 = (expin ^ e) mod q
+ * and then combining the two by the CRT.
+ * 
+ * Two optimisations of this are possible.  First, we can reduce expin
+ * modulo p and q before starting.
+ * 
+ * Second, since we know the factorisation of p and q (trivially derived
+ * from the factorisation of n = p*q), and expin is relatively prime to
+ * both p and q, we can use Euler's theorem, expin^phi(m) = 1 (mod m),
+ * to throw away multiples of phi(p) or phi(q) in e.
+ * Letting ep = e mod phi(p) and
+ *         eq = e mod phi(q)
+ * then combining these two speedups, we only need to evaluate
+ * p2 = ((expin mod p) ^ ep) mod p and
+ * q2 = ((expin mod q) ^ eq) mod q.
+ * 
+ * Now we need to apply the CRT.  Starting with
+ * expout = p2 (mod p) and
+ * expout = q2 (mod q)
+ * we can say that expout = p2 + p * k, and if we assume that 0 <= p2 < p,
+ * then 0 <= expout < p*q for some 0 <= k < q.  Since we want expout = q2
+ * (mod q), then p*k = q2-p2 (mod q).  Since p and q are relatively prime,
+ * p has a multiplicative inverse u mod q.  In other words, u = 1/p (mod q).
+ *
+ * Multiplying by u on both sides gives k = u*(q2-p2) (mod q).
+ * Since we want 0 <= k < q, we can thus find k as
+ * k = (u * (q2-p2)) mod q.
+ * 
+ * Once we have k, evaluating p2 + p * k is easy, and
+ * that gives us the result.
+ * 
+ * In the detailed implementation, there is a temporary, temp, used to
+ * hold intermediate results, p2 is held in expout, and q2 is used as a
+ * temporary in the final step when it is no longer needed.  With that,
+ * you should be able to understand the code below.
+ */
 {
 	unit q2[MAX_UNIT_PRECISION];
 	unit temp[MAX_UNIT_PRECISION];
@@ -1833,9 +1896,9 @@ int mp_modexp_crt(unitptr expout, unitpt
 /*	p2 = [ (expin mod p) ^ ep ] mod p		*/
 	mp_mod(temp,expin,p);		/* temp = expin mod p  */
 	status = mp_modexp(expout,temp,ep,p);
-	if (status < 0)	/* mp_modexp returned an error. */
-	{	mp_init(expout,1);
-		return(status);	/* error return */
+	if (status < 0) {	/* mp_modexp returned an error. */
+		mp_init(expout,1);
+		return status;	/* error return */
 	}
 
 /*	And the same thing for q2 */
@@ -1843,9 +1906,9 @@ int mp_modexp_crt(unitptr expout, unitpt
 /*	q2 = [ (expin mod q) ^ eq ] mod q		*/
 	mp_mod(temp,expin,q);		/* temp = expin mod q  */
 	status = mp_modexp(q2,temp,eq,q);
-	if (status < 0)	/* mp_modexp returned an error. */
-	{	mp_init(expout,1);
-		return(status);	/* error return */
+	if (status < 0) {	/* mp_modexp returned an error. */
+		mp_init(expout,1);
+		return status;	/* error return */
 	}
 
 /*	Now use the multiplicative inverse u to glue together the
@@ -1856,8 +1919,7 @@ int mp_modexp_crt(unitptr expout, unitpt
    but for random results and large q, the odds of (1/q) of it
    being useful do not warrant the test.
 */
-	if (mp_compare(expout,q2) != 0)
-	{
+	if (mp_compare(expout,q2) != 0) {
 #endif
 	/*	Find q2-p2 mod q */
 	if (mp_sub(q2,expout))	/* if the result went negative */
@@ -1874,8 +1936,8 @@ int mp_modexp_crt(unitptr expout, unitpt
 
 	mp_burn(q2);	/* burn the evidence on the stack...*/
 	mp_burn(temp);
-	return(0);	/* normal return */
-}	/* mp_modexp_crt */
+	return 0;	/* normal return */
+} /* mp_modexp_crt */
 
 
 /****************** end of MPI library ****************************/