lbn16.c

包含哈希,对称以及非对称的经典算法 包含经典事例

	xlen = len;
	while (--xlen)
		r = lbnDiv21_16(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) | lbnLshift_16(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 lbnModQ_16
#if defined(BNWORD32) || !BN_SLOW_DIVIDE_32
unsigned
lbnModQ_16(BNWORD16 const *n, unsigned len, unsigned d)
{
	BNWORD16 r;

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

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

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

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

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

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

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

	assert(len > 0);

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

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

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

	while (--len)
		r = lbnDiv21_16(&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 /* lbnModQ_16 */

/*
 * 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.  This is hairy!
 *
 * 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^16) 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_16
BNWORD16
lbnDiv_16(BNWORD16 *q, BNWORD16 *n, unsigned nlen, BNWORD16 *d, unsigned dlen)
{
	BNWORD16 nh,nm,nl;	/* Top three words of the dividend */
	BNWORD16 dh,dl;	/* Top two words of the divisor */
	BNWORD16 qhat;	/* Extimate of quotient word */
	BNWORD16 r;	/* Remainder from quotient estimate division */
	BNWORD16 qhigh;	/* High word of quotient */
	unsigned i;	/* Temp */
	unsigned shift;	/* Bits shifted by normalization */
	unsigned qlen = nlen-dlen; /* Size of quotient (less 1) */
#ifdef mul16_ppmm
	BNWORD16 t16;
#elif defined(BNWORD32)
	BNWORD32 t32;
#else /* use lbnMulN1_16 */
	BNWORD16 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 lbnDiv1_16(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 & ((BNWORD16)1 << 16-1-shift)) == 0) {
			do {
				shift++;
			} while (dh & (BNWORD16)1<<16-1-shift) == 0);
			dh = dh << shift | dl >> (16-shift);
			dl <<= shift;


		}


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

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

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

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

	/* Assert that dh is now normalized */
	assert(dh >> (16-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) {
		assert(nh < dh);
		r = lbnDiv21_16(&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 mul16_ppmm
	mul16_ppmm(nm, t16, qhat, dl);
	if (nm > r || (nm == r && t16 > nl)) {
		/* Decrement qhat and adjust comparison parameters */
		qhat--;
		if ((r += dh) >= dh) {
			nm -= (t16 < dl);
			t16 -= dl;
			if (nm > r || (nm == r && t16 > nl))
				qhat--;
		}
	}
#elif defined(BNWORD32)
	t32 = (BNWORD32)qhat * dl;
	if (t32 > ((BNWORD32)r << 16) + nl) {
		/* Decrement qhat and adjust comparison parameters */
		qhat--;
		if ((r += dh) > dh) {
			t32 -= dl;
			if (t32 > ((BNWORD32)r << 16) + nl)
				qhat--;
		}
	}
#else /* Use lbnMulN1_16 */
	lbnMulN1_16(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
	if (t2high > r || (t2high == r && t2low > nl)) {
		/* Decrement qhat and adjust comparison parameters */
		qhat--;
		if ((r += dh) >= dh) {
			t2high -= (t2low < dl);
			t2low -= dl;
			if (t2high > r || (t2high == r && t2low > nl))
				qhat--;
		}
	}
#endif

	/* Do the multiply and subtract */
	r = lbnMulSub1_16(n, d, dlen, qhat);
	/* If there was a borrow, add back once. */
	if (r > nh) {	/* Borrow? */
		(void)lbnAddN_16(n, d, dlen);
		qhat--;
	}

	/* Remember the first quotient digit. */
	qhigh = qhat;

	/* Now, the main division loop: */
divloop:
	while (qlen--) {

		/* Advance n */
		nh = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));
		BIGLITTLE(++n,--n);
		nm = BIGLITTLE(*(n-dlen),*(n+(dlen-1)));

		if (nh == dh) {
			qhat = ~(BNWORD16)0;
			/* Optimized computation of r = (nh,nm) - qhat * dh */
			r = nh + nm;
			if (r < nh)
				goto subtract;
		} else {
			assert(nh < dh);
			r = lbnDiv21_16(&qhat, nh, nm, dh);
		}

		nl = BIGLITTLE(*(n-(dlen-1)),*(n+(dlen-2)));
#ifdef mul16_ppmm
		mul16_ppmm(nm, t16, qhat, dl);
		if (nm > r || (nm == r && t16 > nl)) {
			/* Decrement qhat and adjust comparison parameters */
			qhat--;
			if ((r += dh) >= dh) {
				nm -= (t16 < dl);
				t16 -= dl;
				if (nm > r || (nm == r && t16 > nl))
					qhat--;
			}
		}
#elif defined(BNWORD32)
		t32 = (BNWORD32)qhat * dl;
		if (t32 > ((BNWORD32)r<<16) + nl) {
			/* Decrement qhat and adjust comparison parameters */
			qhat--;
			if ((r += dh) >= dh) {
				t32 -= dl;
				if (t32 > ((BNWORD32)r << 16) + nl)
					qhat--;
			}
		}
#else /* Use lbnMulN1_16 */
		lbnMulN1_16(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
		if (t2high > r || (t2high == r && t2low > nl)) {
			/* Decrement qhat and adjust comparison parameters */
			qhat--;
			if ((r += dh) >= dh) {
				t2high -= (t2low < dl);
				t2low -= dl;
				if (t2high > r || (t2high == r && t2low > nl))
					qhat--;
			}
		}
#endif

		/*
		 * As a point of interest, note that it is not worth checking
		 * for qhat of 0 or 1 and installing special-case code.  These
		 * occur with probability 2^-16, so spending 1 cycle to check
		 * for them is only worth it if we save more than 2^15 cycles,
		 * and a multiply-and-subtract for numbers in the 1024-bit
		 * range just doesn't take that long.
		 */
subtract:
		/*
		 * n points to the least significant end of the substring
		 * of n to be subtracted from.  qhat is either exact or
		 * one too large.  If the subtract gets a borrow, it was
		 * one too large and the divisor is added back in.  It's
		 * a dlen+1 word add which is guaranteed to produce a
		 * carry out, so it can be done very simply.
		 */
		r = lbnMulSub1_16(n, d, dlen, qhat);
		if (r > nh) {	/* Borrow? */
			(void)lbnAddN_16(n, d, dlen);
			qhat--;
		}
		/* Store the quotient digit */
		BIGLITTLE(*q++,*--q) = qhat;
	}
	/* Tah dah! */

	if (shift) {
		lbnRshift_16(d, dlen, shift);
		lbnRshift_16(n, dlen, shift);
	}

	return qhigh;