bni16.c

PGP—Pretty Good Privacy

	BNWORD16 *pp;
	BNWORD16 t;
	unsigned carry;
	unsigned i, j;

	/* Special case of zero */
	if (!len)
		return;

	/*
	 * This computes directly into the high half of prod, so just
	 * shift the pointer and consider prod only "len" elements long
	 * for the rest of the code.
	 */
	BIGLITTLE(prod -= len, prod += len);

	/* Pass 1 - compute Montgomery multipliers */
	/* First iteration can have certain simplifications. */
	x = (BNWORD32)BIGLITTLE(num1[-1] * num2[-1], num1[0] * num2[0]);
	BIGLITTLE(prod[-1], prod[0]) = t = inv * (BNWORD16)x;
	y = (BNWORD32)t * BIGLITTLE(mod[-1],mod[0]);
	x += y;
	/* Note: GCC 2.6.3 has a bug if you try to eliminate "carry" */
	carry = (x < y);
	pgpAssert((BNWORD16)x == 0);
	x = x >> 16 | (BNWORD32)carry << 16;

	for (i = 1; i < len; i++) {
		carry = 0;
		p1 = num1;
		p2 = BIGLITTLE(num2-i-1,num2+i+1);
		pp = prod;
		pm = BIGLITTLE(mod-i-1,mod+i+1);
		for (j = 0; j < i; j++) {
			y = (BNWORD32)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
			x += y;
			carry += (x < y);
			y = (BNWORD32)BIGLITTLE(*--pp * *pm++, *pp++ * *--pm);
			x += y;
			carry += (x < y);
		}
		y = (BNWORD32)BIGLITTLE(p1[-1] * p2[0], p1[0] * p2[-1]);
		x += y;
		carry += (x < y);
		pgpAssert(BIGLITTLE(pp == prod-i, pp == prod+i));
		BIGLITTLE(pp[-1], pp[0]) = t = inv * (BNWORD16)x;
		pgpAssert(BIGLITTLE(pm == mod-1, pm == mod+1));
		y = (BNWORD32)t * BIGLITTLE(pm[0],pm[-1]);
		x += y;
		carry += (x < y);
		pgpAssert((BNWORD16)x == 0);
		x = x >> 16 | (BNWORD32)carry << 16;
	}

	/* Pass 2 - compute reduced product and store */
	for (i = 1; i < len; i++) {
		carry = 0;
		p1 = BIGLITTLE(num1-i,num1+i);
		p2 = BIGLITTLE(num2-len,num2+len);
		pm = BIGLITTLE(mod-i,mod+i);
		pp = BIGLITTLE(prod-len,prod+len);
		for (j = i; j < len; j++) {
			y = (BNWORD32)BIGLITTLE(*--p1 * *p2++, *p1++ * *--p2);
			x += y;
			carry += (x < y);
			y = (BNWORD32)BIGLITTLE(*--pm * *pp++, *pm++ * *--pp);
			x += y;
			carry += (x < y);
		}
		pgpAssert(BIGLITTLE(pm == mod-len, pm == mod+len));
		pgpAssert(BIGLITTLE(pp == prod-i, pp == prod+i));
		BIGLITTLE(pp[0],pp[-1]) = (BNWORD16)x;
		x = (x >> 16) | (BNWORD32)carry << 16;
	}

	/* Last round of second half, simplified. */
	BIGLITTLE(*(prod-len),*(prod+len-1)) = (BNWORD16)x;
	carry = (x >> 16);

	while (carry)
		carry -= bniSubN_16(prod, mod, len);
	while (bniCmp_16(prod, mod, len) >= 0)
		(void)bniSubN_16(prod, mod, len);
}
/* Suppress later definition */
#define bniMontMul_16 bniMontMul_16
#endif

#if !defined(bniSquare_16) && defined(BNWORD32) && PRODUCT_SCAN
/*
 * Trial code for product-scanning squaring.  This seems to slow the C
 * code down rather than speed it up.
 */
void
bniSquare_16(BNWORD16 *prod, BNWORD16 const *num, unsigned len)
{
	BNWORD32 x, y, z;
	BNWORD16 const *p1, *p2;
	unsigned carry;
	unsigned i, j;

	/* Special case of zero */
	if (!len)
		return;

	/* Word 0 of product */
	x = (BNWORD32)BIGLITTLE(num[-1] * num[-1], num[0] * num[0]);
	BIGLITTLE(*--prod, *prod++) = (BNWORD16)x;
	x >>= 16;

	/* Words 1 through len-1 */
	for (i = 1; i < len; i++) {
		carry = 0;
		y = 0;
		p1 = num;
		p2 = BIGLITTLE(num-i-1,num+i+1);
		for (j = 0; j < (i+1)/2; j++) {
			BIG(z = (BNWORD32)*--p1 * *p2++;)
			LITTLE(z = (BNWORD32)*p1++ * *--p2;)
			y += z;
			carry += (y < z);
		}
		y += z = y;
		carry += carry + (y < z);
		if ((i & 1) == 0) {
			pgpAssert(BIGLITTLE(p1-1 == p2, p1 == p2-1));
			BIG(z = (BNWORD32)*p2 * *p2;)
			LITTLE(z = (BNWORD32)*p1 * *p1;)
			y += z;
			carry += (y < z);
		}
		x += y;
		carry += (x < y);
		BIGLITTLE(*--prod,*prod++) = (BNWORD16)x;
		x = (x >> 16) | (BNWORD32)carry << 16;
	}
	/* Words len through 2*len-2 */
	for (i = 1; i < len; i++) {
		carry = 0;
		y = 0;
		p1 = BIGLITTLE(num-i,num+i);
		p2 = BIGLITTLE(num-len,num+len);
		for (j = 0; j < (len-i)/2; j++) {
			BIG(z = (BNWORD32)*--p1 * *p2++;)
			LITTLE(z = (BNWORD32)*p1++ * *--p2;)
			y += z;
			carry += (y < z);
		}
		y += z = y;
		carry += carry + (y < z);
		if ((len-i) & 1) {
			pgpAssert(BIGLITTLE(p1-1 == p2, p1 == p2-1));
			BIG(z = (BNWORD32)*p2 * *p2;)
			LITTLE(z = (BNWORD32)*p1 * *p1;)
			y += z;
			carry += (y < z);
		}
		x += y;
		carry += (x < y);
		BIGLITTLE(*--prod,*prod++) = (BNWORD16)x;
		x = (x >> 16) | (BNWORD32)carry << 16;
	}
	
	/* Word 2*len-1 */
	BIGLITTLE(*--prod,*prod) = (BNWORD16)x;
}
/* Suppress later definition */
#define bniSquare_16 bniSquare_16
#endif

/*
 * Square a number, using optimized squaring to reduce the number of
 * primitive multiples that are executed.  There may not be any
 * overlap of the input and output.
 *
 * Technique: Consider the partial products in the multiplication
 * of "abcde" by itself:
 *
 *               a  b  c  d  e
 *            *  a  b  c  d  e
 *          ==================
 *              ae be ce de ee
 *           ad bd cd dd de
 *        ac bc cc cd ce
 *     ab bb bc bd be
 *  aa ab ac ad ae
 *
 * Note that everything above the main diagonal:
 *              ae be ce de = (abcd) * e
 *           ad bd cd       = (abc) * d
 *        ac bc             = (ab) * c
 *     ab                   = (a) * b
 *
 * is a copy of everything below the main diagonal:
 *                       de
 *                 cd ce
 *           bc bd be
 *     ab ac ad ae
 *
 * Thus, the sum is 2 * (off the diagonal) + diagonal.
 *
 * This is accumulated beginning with the diagonal (which
 * consist of the squares of the digits of the input), which is then
 * divided by two, the off-diagonal added, and multiplied by two
 * again.  The low bit is simply a copy of the low bit of the
 * input, so it doesn't need special care.
 *
 * TODO: Merge the shift by 1 with the squaring loop.
 * TODO: Use Karatsuba.  (a*W+b)^2 = a^2 * (W^2+W) + b^2 * (W+1) - (a-b)^2 * W.
 */
#ifndef bniSquare_16
void
bniSquare_16(BNWORD16 *prod, BNWORD16 const *num, unsigned len)
{
	BNWORD16 t;
	BNWORD16 *prodx = prod;		/* Working copy of the argument */
	BNWORD16 const *numx = num;	/* Working copy of the argument */
	unsigned lenx = len;		/* Working copy of the argument */

	if (!len)
		return;

	/* First, store all the squares */
	while (lenx--) {
#ifdef mul16_ppmm
		BNWORD16 ph, pl;
		t = BIGLITTLE(*--numx,*numx++);
		mul16_ppmm(ph,pl,t,t);
		BIGLITTLE(*--prodx,*prodx++) = pl;
		BIGLITTLE(*--prodx,*prodx++) = ph;
#elif defined(BNWORD32) /* use BNWORD32 */
		BNWORD32 p;
		t = BIGLITTLE(*--numx,*numx++);
		p = (BNWORD32)t * t;
		BIGLITTLE(*--prodx,*prodx++) = (BNWORD16)p;
		BIGLITTLE(*--prodx,*prodx++) = (BNWORD16)(p>>16);
#else	/* Use bniMulN1_16 */
		t = BIGLITTLE(numx[-1],*numx);
		bniMulN1_16(prodx, numx, 1, t);
		BIGLITTLE(--numx,numx++);
		BIGLITTLE(prodx -= 2, prodx += 2);
#endif
	}
	/* Then, shift right 1 bit */
	(void)bniRshift_16(prod, 2*len, 1);

	/* Then, add in the off-diagonal sums */
	lenx = len;
	numx = num;
	prodx = prod;
	while (--lenx) {
		t = BIGLITTLE(*--numx,*numx++);
		BIGLITTLE(--prodx,prodx++);
		t = bniMulAdd1_16(prodx, numx, lenx, t);
		bniAdd1_16(BIGLITTLE(prodx-lenx,prodx+lenx), lenx+1, t);
		BIGLITTLE(--prodx,prodx++);
	}

	/* Shift it back up */
	bniDouble_16(prod, 2*len);

	/* And set the low bit appropriately */
	BIGLITTLE(prod[-1],prod[0]) |= BIGLITTLE(num[-1],num[0]) & 1;
}
#endif /* !bniSquare_16 */

/*
 * bniNorm_16 - given a number, return a modified length such that the
 * most significant digit is non-zero.  Zero-length input is okay.
 */
#ifndef bniNorm_16
unsigned
bniNorm_16(BNWORD16 const *num, unsigned len)
{
	BIGLITTLE(num -= len,num += len);
	while (len && BIGLITTLE(*num++,*--num) == 0)
		--len;
	return len;
}
#endif /* bniNorm_16 */

/*
 * bniBits_16 - return the number of significant bits in the array.
 * It starts by normalizing the array.  Zero-length input is okay.
 * Then assuming there's anything to it, it fetches the high word,
 * generates a bit length by multiplying the word length by 16, and
 * subtracts off 16/2, 16/4, 16/8, ... bits if the high bits are clear.
 */
#ifndef bniBits_16
unsigned
bniBits_16(BNWORD16 const *num, unsigned len)
{
	BNWORD16 t;
	unsigned i;

	len = bniNorm_16(num, len);
	if (len) {
		t = BIGLITTLE(*(num-len),*(num+(len-1)));
		pgpAssert(t);
		len *= 16;
		i = 16/2;
		do {
			if (t >> i)
				t >>= i;
			else
				len -= i;
		} while ((i /= 2) != 0);
	}
	return len;
}
#endif /* bniBits_16 */

/*
 * 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 16 bits, which a general 32-bit divide
 * in a compiler's run-time library can't do.
 */
#ifndef BN_SLOW_DIVIDE_32
/* Assume that divisors of more than thirty-two bits are slow */
#define BN_SLOW_DIVIDE_32 (32 > 0x20)
#endif

/*
 * Return (nh<<16|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_16
#if defined(BNWORD32) && !BN_SLOW_DIVIDE_32
BNWORD16
bniDiv21_16(BNWORD16 *q, BNWORD16 nh, BNWORD16 nl, BNWORD16 d)
{
	BNWORD32 n = (BNWORD32)nh << 16 | nl;

	/* Divisor must be normalized */
	pgpAssert(d >> (16-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 16/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) >> 16/2 )
#define lowhalf(x) ( (x) & (((BNWORD16)1 << 16/2)-1) )
BNWORD16
bniDiv21_16(BNWORD16 *q, BNWORD16 nh, BNWORD16 nl, BNWORD16 d)
{
	BNWORD16 dh = highhalf(d), dl = lowhalf(d);
	BNWORD16 qh, ql, prod, r;

	/* Divisor must be normalized */
	pgpAssert((d >> (16-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 << (16/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 << (16/2)) | lowhalf(nl);
	if (r < prod) {
		--ql; r += d;
		if (r >= d && r < prod) {
			--ql; r += d;
		}
	}
	r -= prod;