bignumber.cpp

绝对好东东

{
	BNWORD32 t;
	unsigned i;

	len = bniNorm_32(num, len);
	if (len) {
		t = BIGLITTLE(*(num-len),*(num+(len-1)));
		ASSERT(t);
		len *= 32;
		i = 32/2;
		do {
			if (t >> i)
				t >>= i;
			else
				len -= i;
		} while ((i /= 2) != 0);
	}
	return len;
}
#endif /* bniBits_32 */

/*
 * If defined, use hand-rolled divide rather than compiler's native.
 * If the machine doesn't do it in line, the manual code is probably
 * faster, since it can assume normalization and the fact that the
 * quotient will fit into 32 bits, which a general 64-bit divide
 * in a compiler's run-time library can't do.
 */
#ifndef BN_SLOW_DIVIDE_64
/* Assume that divisors of more than thirty-two bits are slow */
#define BN_SLOW_DIVIDE_64 (64 > 0x20)
#endif

/*
 * Return (nh<<32|nl) % d, and place the quotient digit into *q.
 * It is guaranteed that nh < d, and that d is normalized (with its high
 * bit set).  If we have a double-width type, it's easy.  If not, ooh,
 * yuk!
 */
#ifndef bniDiv21_32
#if defined(BNWORD64) && !BN_SLOW_DIVIDE_64
BNWORD32
bniDiv21_32(BNWORD32 *q, BNWORD32 nh, BNWORD32 nl, BNWORD32 d)
{
	BNWORD64 n = (BNWORD64)nh << 32 | nl;

	/* Divisor must be normalized */
	ASSERT(d >> (32-1) == 1);

	*q = n / d;
	return n % d;
}
#else
/*
 * This is where it gets ugly.
 *
 * Do the division in two halves, using Algorithm D from section 4.3.1
 * of Knuth.  Note Theorem B from that section, that the quotient estimate
 * is never more than the true quotient, and is never more than two
 * too low.
 *
 * The mapping onto conventional long division is (everything a half word):
 *        _____________qh___ql_
 * dh dl ) nh.h nh.l nl.h nl.l
 *             - (qh * d)
 *            -----------
 *              rrrr rrrr nl.l
 *                  - (ql * d)
 *                -----------
 *                  rrrr rrrr
 *
 * The implicit 3/2-digit d*qh and d*ql subtractors are computed this way:
 *   First, estimate a q digit so that nh/dh works.  Subtracting qh*dh from
 *   the (nh.h nh.l) list leaves a 1/2-word remainder r.  Then compute the
 *   low part of the subtractor, qh * dl.   This also needs to be subtracted
 *   from (nh.h nh.l nl.h) to get the final remainder.  So we take the
 *   remainder, which is (nh.h nh.l) - qh*dl, shift it and add in nl.h, and
 *   try to subtract qh * dl from that.  Since the remainder is 1/2-word
 *   long, shifting and adding nl.h results in a single word r.
 *   It is possible that the remainder we're working with, r, is less than
 *   the product qh * dl, if we estimated qh too high.  The estimation
 *   technique can produce a qh that is too large (never too small), leading
 *   to r which is too small.  In that case, decrement the digit qh, add
 *   shifted dh to r (to correct for that error), and subtract dl from the
 *   product we're comparing r with.  That's the "correct" way to do it, but
 *   just adding dl to r instead of subtracting it from the product is
 *   equivalent and a lot simpler.  You just have to watch out for overflow.
 *
 *   The process is repeated with (rrrr rrrr nl.l) for the low digit of the
 *   quotient ql.
 *
 * The various uses of 32/2 for shifts are because of the note about
 * automatic editing of this file at the very top of the file.
 */
#define highhalf(x) ( (x) >> 32/2 )
#define lowhalf(x) ( (x) & (((BNWORD32)1 << 32/2)-1) )
BNWORD32
bniDiv21_32(BNWORD32 *q, BNWORD32 nh, BNWORD32 nl, BNWORD32 d)
{
	BNWORD32 dh = highhalf(d), dl = lowhalf(d);
	BNWORD32 qh, ql, prod, r;

	/* Divisor must be normalized */
	ASSERT((d >> (32-1)) == 1);

	/* Do first half-word of division */
	qh = nh / dh;
	r = nh % dh;
	prod = qh * dl;

	/*
	 * Add next half-word of numerator to remainder and correct.
	 * qh may be up to two too large.
	 */
	r = (r << (32/2)) | highhalf(nl);
	if (r < prod) {
		--qh; r += d;
		if (r >= d && r < prod) {
			--qh; r += d; 
		}
	}
	r -= prod;

	/* Do second half-word of division */
	ql = r / dh;
	r = r % dh;
	prod = ql * dl;

	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;

	ASSERT(len > 0);
	ASSERT(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);
	ASSERT((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;

	ASSERT(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;

	ASSERT(len > 0);

	shift = 0;
	r = d;
	i = 32;
	while (i /= 2) {
		if (r >> i)
			r >>= i;
		else
			shift += i;
	}
	ASSERT(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

	ASSERT(dlen);
	ASSERT(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);