bni16.c

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

* times the Montgomery multipliers. The results of this multiply
* are stored.
*/
static void
bniMontMul_16(BNWORD16 *prod, BNWORD16 const *num1, BNWORD16 const *num2,
			BNWORD16 const *mod, unsigned len, BNWORD16 inv)
	{
			BNWORD32 x, y;
			BNWORD16 const *p1, *p2, *pm;
			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;