bni32.c

著名的加密软件的应用于电子邮件中

r = (r << (32/2)) | lowhalf(nl);
if (r < prod) {
			--ql; r += d;
			if (r >= d && r < prod) {
				 --ql; r += d;
			}
}
r -= prod;

*q = (qh << (32/2)) | ql;

return r;
}
#endif
#endif /* bniDiv21_32 */


/*
* In the division functions, the dividend and divisor are referred to
* as "n" and "d", which stand for "numerator" and "denominator".
*
* The quotient is (nlen-dlen+1) digits long. It may be overlapped with
* the high (nlen-dlen) words of the dividend, but one extra word is needed
* on top to hold the top word.
*/

/*
* Divide an n-word number by a 1-word number, storing the remainder
* and n-1 words of the n-word quotient. The high word is returned.
* It IS legal for rem to point to the same address as n, and for
* q to point one word higher.
*
* TODO: If BN_SLOW_DIVIDE_64, add a divnhalf_32 which uses 32-bit
* dividends if the divisor is half that long.
* TODO: Shift the dividend on the fly to avoid the last division and
* instead have a remainder that needs shifting.
* TODO: Use reciprocals rather than dividing.
	*/
#ifndef bniDiv1_32
BNWORD32
bniDiv1_32(BNWORD32 *q, BNWORD32 *rem, BNWORD32 const *n, unsigned len,
			BNWORD32 d)
	{
			unsigned shift;
			unsigned xlen;
			BNWORD32 r;
			BNWORD32 qhigh;

pgpAssert(len > 0);
pgpAssert(d);

if (len == 1) {
r = *n;
*rem = r%d;
return r/d;
}

			shift = 0;
			r = d;
			xlen = 32/2;
			do {
					if (r >> xlen)
						r >>= xlen;
					else
						shift += xlen;
			} while ((xlen /= 2) != 0);
			pgpAssert((d >> (32-1-shift)) == 1);
			d <<= shift;

			BIGLITTLE(q -= len-1,q += len-1);
			BIGLITTLE(n -= len,n += len);

			r = BIGLITTLE(*n++,*--n);
			if (r < d) {
				 qhigh = 0;
			} else {
			qhigh = r/d;
				 r %= d;
			}

			xlen = len;
			while (--xlen)
				 r = bniDiv21_32(BIGLITTLE(q++,--q), r,
				 		 BIGLITTLE(*n++,*--n), d);

			/*
			* Final correction for shift - shift the quotient up "shift"
			* bits, and merge in the extra bits of quotient. Then reduce
			* the final remainder mod the real d.
			*/
			if (shift) {
					d >>= shift;
					qhigh = (qhigh << shift) | bniLshift_32(q, len-1, shift);
					BIGLITTLE(q[-1],*q) |= r/d;
					r %= d;
			}
			*rem = r;

			return qhigh;
	}
	#endif

/*
* This function performs a "quick" modulus of a number with a divisor
* d which is guaranteed to be at most sixteen bits, i.e. less than 65536.
* This applies regardless of the word size the library is compiled with.
*
* This function is important to prime generation, for sieving.
*/
#ifndef bniModQ_32
/* If there's a custom bniMod21_32, no normalization needed */
#ifdef bniMod21_32
unsigned
bniModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
	{
			unsigned i, shift;
			BNWORD32 r;

			pgpAssert(len > 0);

			BIGLITTLE(n -= len,n += len);

		/* Try using a compare to avoid the first divide */
			r = BIGLITTLE(*n++,*--n);
			if (r >= d)
				 r %= d;
			while (--len)
				 r = bniMod21_32(r, BIGLITTLE(*n++,*--n), d);

			return r;
	}
#elif defined(BNWORD64) && !BN_SLOW_DIVIDE_64
unsigned
bniModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
{
			BNWORD32 r;

if (!--len)
	return BIGLITTLE(n[-1],n[0]) % d;

BIGLITTLE(n -= len,n += len);
r = BIGLITTLE(n[-1],n[0]);

do {
	r = (BNWORD32)((((BNWORD64)r<<32) | BIGLITTLE(*n++,*--n)) % d);
} while (--len);

return r;
}
#elif 32 >= 0x20
/*
* If the single word size can hold 65535*65536, then this function
* is avilable.
*/
#ifndef highhalf
#define highhalf(x) ( (x) >> 32/2 )
#define lowhalf(x) ( (x) & ((1 << 32/2)-1) )
#endif
unsigned
bniModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
	{
			BNWORD32 r, x;

			BIGLITTLE(n -= len,n += len);

			r = BIGLITTLE(*n++,*--n);
			while (--len) {
				 x = BIGLITTLE(*n++,*--n);
				 r = (r%d << 32/2) | highhalf(x);
				 r = (r%d << 32/2) | lowhalf(x);
			}

			return r%d;
	}
#else
/* Default case - use bniDiv21_32 */
unsigned
bniModQ_32(BNWORD32 const *n, unsigned len, unsigned d)
	{
			unsigned i, shift;
			BNWORD32 r;
			BNWORD32 q;

			pgpAssert(len > 0);

			shift = 0;
			r = d;
			i = 32;
			while (i /= 2) {
					if (r >> i)
						r >>= i;
					else
						shift += i;
			}
			pgpAssert(d >> (32-1-shift) == 1);
			d <<= shift;

			BIGLITTLE(n -= len,n += len);

			r = BIGLITTLE(*n++,*--n);
			if (r >= d)
				 r %= d;

			while (--len)
				 r = bniDiv21_32(&q, r, BIGLITTLE(*n++,*--n), d);

		/*
		* Final correction for shift - shift the quotient up "shift"
		* bits, and merge in the extra bits of quotient. Then reduce
		* the final remainder mod the real d.
		*/
		if (shift)
			r %= d >> shift;

			return r;
	}
	#endif
#endif /* bniModQ_32 */

/*
* Reduce n mod d and return the quotient. That is, find:
* q = n / d;
* n = n % d;
* d is altered during the execution of this subroutine by normalizing it.
* It must already have its most significant word non-zero; it is shifted
* so its most significant bit is non-zero.
*
* The quotient q is nlen-dlen+1 words long. To make it possible to
* overlap the quptient with the input (you can store it in the high dlen
* words), the high word of the quotient is *not* stored, but is returned.
* (If all you want is the remainder, you don't care about it, anyway.)
*
* This uses algorithm D from Knuth (4.3.1), except that we do binary
* (shift) normalization of the divisor. WARNING: This is hairy!
*
* This function is used for some modular reduction, but it is not used in
* the modular exponentiation loops; they use Montgomery form and the
* corresponding, more efficient, Montgomery reduction. This code
* is needed for the conversion to Montgomery form, however, so it
* has to be here and it might as well be reasonably efficient.
*
* The overall operation is as follows ("top" and "up" refer to the
* most significant end of the number; "bottom" and "down", the least):
*
* - Shift the divisor up until the most significant bit is set.
* - Shift the dividend up the same amount. This will produce the
* correct quotient, and the remainder can be recovered by shifting
* it back down the same number of bits. This may produce an overflow
* word, but the word is always strictly less than the most significant
* divisor word.
* - Estimate the first quotient digit qhat:
* - First take the top two words (one of which is the overflow) of the
* dividend and divide by the top word of the divisor:
* qhat = (nh,nm)/dh. This qhat is >= the correct quotient digit
* and, since dh is normalized, it is at most two over.
* - Second, correct by comparing the top three words. If
*	(dh,dl) * qhat > (nh,nm,ml), decrease qhat and try again.
*	The second iteration can be simpler because there can't be a third.
*	The computation can be simplified by subtracting dh*qhat from
*	both sides, suitably shifted. This reduces the left side to
*	dl*qhat. On the right, (nh,nm)-dh*qhat is simply the
*	remainder r from (nh,nm)%dh, so the right is (r,nl).
*	This produces qhat that is almost always correct and at
*	most (prob ~ 2/2^32) one too high.
* - Subtract qhat times the divisor (suitably shifted) from the dividend.
* If there is a borrow, qhat was wrong, so decrement it
* and add the divisor back in (once).
* - Store the final quotient digit qhat in the quotient array q.
*
* Repeat the quotient digit computation for successive digits of the
* quotient until the whole quotient has been computed. Then shift the
* divisor and the remainder down to correct for the normalization.
*
* TODO: Special case 2-word divisors.
* TODO: Use reciprocals rather than dividing.
	*/
#ifndef divn_32
BNWORD32
bniDiv_32(BNWORD32 *q, BNWORD32 *n, unsigned nlen, BNWORD32 *d, unsigned dlen)
{
BNWORD32 nh,nm,nl;	/* Top three words of the dividend */
BNWORD32 dh,dl;	/* Top two words of the divisor */
BNWORD32 qhat;	/* Extimate of quotient word */
BNWORD32 r;	 /* Remainder from quotient estimate division */
	BNWORD32 qhigh;	/* High word of quotient */
	unsigned i;	/* Temp */
	unsigned shift;	/* Bits shifted by normalization */
	unsigned qlen = nlen-dlen; /* Size of quotient (less 1) */
#ifdef mul32_ppmm
	BNWORD32 t32;
#elif defined(BNWORD64)
	BNWORD64 t64;
#else /* use bniMulN1_32 */
	BNWORD32 t2[2];
#define t2high BIGLITTLE(t2[0],t2[1])
#define t2low BIGLITTLE(t2[1],t2[0])
	#endif

		pgpAssert(dlen);
		pgpAssert(nlen >= dlen);

		/*
			* Special cases for short divisors. The general case uses the
			* top top 2 digits of the divisor (d) to estimate a quotient digit,
			* so it breaks if there are fewer digits available. Thus, we need
			* special cases for a divisor of length 1. A divisor of length
			* 2 can have a *lot* of administrivia overhead removed removed,
			* so it's probably worth special-casing that case, too.
			*/
			if (dlen == 1)
				 return bniDiv1_32(q, BIGLITTLE(n-1,n), n, nlen,
						BIGLITTLE(d[-1],d[0]));

#if 0
/*
* @@@ This is not yet written... The general loop will do,
* albeit less efficiently
*/
if (dlen == 2) {
				/*
			* divisor two digits long:
			* use the 3/2 technique from Knuth, but we know
			* it's exact.
			*/
			dh = BIGLITTLE(d[-1],d[0]);
			dl = BIGLITTLE(d[-2],d[1]);
			shift = 0;
			if ((sh & ((BNWORD32)1 << 32-1-shift)) == 0) {
						do {
							shift++;
						} while (dh & (BNWORD32)1<<32-1-shift) == 0);
						dh = dh << shift | dl >> (32-shift);
						dl <<= shift;


}


				 for (shift = 0; (dh & (BNWORD32)1 << 32-1-shift)) == 0;
									shift++)
				 	;
				 if (shift) {
}
dh = dh << shift | dl >> (32-shift);
shift = 0;
while (dh
}
#endif

			dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
			pgpAssert(dh);

			/* Normalize the divisor */
			shift = 0;
			r = dh;
			i = 32/2;
			do {
					if (r >> i)
						r >>= i;
					else
						shift += i;
			} while ((i /= 2) != 0);

			nh = 0;
			if (shift) {
				 bniLshift_32(d, dlen, shift);
				 dh = BIGLITTLE(*(d-dlen),*(d+(dlen-1)));
				 nh = bniLshift_32(n, nlen, shift);
			}

			/* Assert that dh is now normalized */
			pgpAssert(dh >> (32-1));

			/* Also get the second-most significant word of the divisor */
			dl = BIGLITTLE(*(d-(dlen-1)),*(d+(dlen-2)));

			/*
			* Adjust pointers: n to point to least significant end of first
			* first subtract, and q to one the most-significant end of the
			* quotient array.
			*/
			BIGLITTLE(n -= qlen,n += qlen);
			BIGLITTLE(q -= qlen,q += qlen);

			/* Fetch the most significant stored word of the dividend */
			nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));

			/*
			* Compute the first digit of the quotient, based on the
			* first two words of the dividend (the most significant of which
			* is the overflow word h).
			*/
			if (nh) {
				 pgpAssert(nh < dh);
				 r = bniDiv21_32(&qhat, nh, nm, dh);
			} else if (nm >= dh) {
				 qhat = nm/dh;
				 r = nm % dh;
			} else {	 /* Quotient is zero */
				 qhigh = 0;
	goto divloop;
}

/* Now get the third most significant word of the dividend */
nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));

			/*
			* Correct qhat, the estimate of quotient digit.
			* qhat can only be high, and at most two words high,
			* so the loop can be unrolled and abbreviated.
			*/
#ifdef mul32_ppmm
			mul32_ppmm(nm, t32, qhat, dl);
			if (nm > r || (nm == r && t32 > nl)) {
					/* Decrement qhat and adjust comparison parameters */
					qhat--;
					if ((r += dh) >= dh) {
							nm -= (t32 < dl);
							t32 -= dl;
							if (nm > r || (nm == r && t32 > nl))
								 qhat--;
					}
			}
#elif defined(BNWORD64)
		t64 = (BNWORD64)qhat * dl;
			if (t64 > ((BNWORD64)r << 32) + nl) {