bni16.c

PGP—Pretty Good Privacy

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

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


/*
 * 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_32, add a divnhalf_16 which uses 16-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_16
BNWORD16
bniDiv1_16(BNWORD16 *q, BNWORD16 *rem, BNWORD16 const *n, unsigned len,
	BNWORD16 d)
{
	unsigned shift;
	unsigned xlen;
	BNWORD16 r;
	BNWORD16 qhigh;

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

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

	shift = 0;
	r = d;
	xlen = 16/2;
	do {
		if (r >> xlen)
			r >>= xlen;
		else
			shift += xlen;
	} while ((xlen /= 2) != 0);
	pgpAssert((d >> (16-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_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) | bniLshift_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 bniModQ_16
/* If there's a custom bniMod21_16, no normalization needed */
#ifdef bniMod21_16
unsigned
bniModQ_16(BNWORD16 const *n, unsigned len, unsigned d)
{
	unsigned i, shift;
	BNWORD16 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_16(r, BIGLITTLE(*n++,*--n), d);

	return r;
}
#elif defined(BNWORD32) && !BN_SLOW_DIVIDE_32
unsigned
bniModQ_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
bniModQ_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 bniDiv21_16 */
unsigned
bniModQ_16(BNWORD16 const *n, unsigned len, unsigned d)
{
	unsigned i, shift;
	BNWORD16 r;
	BNWORD16 q;

	pgpAssert(len > 0);

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

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

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

	while (--len)
		r = bniDiv21_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 /* bniModQ_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.  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^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
bniDiv_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 bniMulN1_16 */
	BNWORD16 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_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)));
	pgpAssert(dh);

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

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

	/* Assert that dh is now normalized */
	pgpAssert(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) {
		pgpAssert(nh < dh);
		r = bniDiv21_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 bniMulN1_16 */
	bniMulN1_16(BIGLITTLE(t2+2,t2), &dl, 1, qhat);
	if (t2high > r || (t2high == r && t2low > nl)) {