mpi.c

matrix ssl代码

		return;
	}

	{
		register mp_digit *bottom, *top;

/*
		Shift the digits down
 */
		/* bottom */
		bottom = a->dp;

		/* top [offset into digits] */
		top = a->dp + b;

/*
		This is implemented as a sliding window where the window is b-digits long
		and digits from the top of the window are copied to the bottom.
		
		 e.g.

		b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->
		/\                   |      ---->
		\-------------------/      ---->
 */
		for (x = 0; x < (a->used - b); x++) {
			*bottom++ = *top++;
		}

/*
		Zero the top digits
 */
		for (; x < a->used; x++) {
			*bottom++ = 0;
		}
	}

/*
	Remove excess digits
 */
	a->used -= b;
}

/******************************************************************************/
/* 
	Low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9
 */
int32 s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
{
	int32		olduse, res, min, max;

/*
	Find sizes
 */
	min = b->used;
	max = a->used;

/*
	init result
 */
	if (c->alloc < max) {
		if ((res = mp_grow (c, max)) != MP_OKAY) {
			return res;
		}
	}
	olduse = c->used;
	c->used = max;

	{
		register mp_digit u, *tmpa, *tmpb, *tmpc;
		register int32 i;

/*
		alias for digit pointers
 */
		tmpa = a->dp;
		tmpb = b->dp;
		tmpc = c->dp;

/*
		set carry to zero
 */
		u = 0;
		for (i = 0; i < min; i++) {
			/* T[i] = A[i] - B[i] - U */
			*tmpc = *tmpa++ - *tmpb++ - u;

/*
			U = carry bit of T[i]
			Note this saves performing an AND operation since if a carry does occur it
			will propagate all the way to the MSB.  As a result a single shift
			is enough to get the carry
 */
			u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));

			/* Clear carry from T[i] */
			*tmpc++ &= MP_MASK;
		}

/*
		Now copy higher words if any, e.g. if A has more digits than B
 */
		for (; i < max; i++) {
			/* T[i] = A[i] - U */
			*tmpc = *tmpa++ - u;

			/* U = carry bit of T[i] */
			u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));

			/* Clear carry from T[i] */
			*tmpc++ &= MP_MASK;
		}

/*
		Clear digits above used (since we may not have grown result above)
 */
		for (i = c->used; i < olduse; i++) {
			*tmpc++ = 0;
		}
	}

	mp_clamp (c);
	return MP_OKAY;
}
/******************************************************************************/
/*
	integer signed division. 

	c*b + d == a [e.g. a/b, c=quotient, d=remainder]
	HAC pp.598 Algorithm 14.20

	Note that the description in HAC is horribly incomplete.  For example,
	it doesn't consider the case where digits are removed from 'x' in the inner
	loop.  It also doesn't consider the case that y has fewer than three
	digits, etc..

	The overall algorithm is as described as 14.20 from HAC but fixed to
	treat these cases.
 */
#ifdef MP_DIV_SMALL
int32 mp_div(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c, mp_int * d)
{
	mp_int	ta, tb, tq, q;
	int32		res, n, n2;

/*
	is divisor zero ?
 */
	if (mp_iszero (b) == 1) {
		return MP_VAL;
	}

/*
	if a < b then q=0, r = a
 */
	if (mp_cmp_mag (a, b) == MP_LT) {
		if (d != NULL) {
			res = mp_copy (a, d);
		} else {
			res = MP_OKAY;
		}
		if (c != NULL) {
			mp_zero (c);
		}
		return res;
	}
	
/*
	init our temps
 */
	if ((res = _mp_init_multi(pool, &ta, &tb, &tq, &q, NULL, NULL, NULL, NULL) != MP_OKAY)) {
		return res;
	}

/*
	tq = 2^n,  tb == b*2^n
 */
	mp_set(&tq, 1);
	n = mp_count_bits(a) - mp_count_bits(b);
	if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
			((res = mp_abs(b, &tb)) != MP_OKAY) || 
			((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
			((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
      goto __ERR;
	}
/* old
	if (((res = mp_copy(a, &ta)) != MP_OKAY) ||
		((res = mp_copy(b, &tb)) != MP_OKAY) || 
		((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
		((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
			goto LBL_ERR;
	}
*/
	while (n-- >= 0) {
		if (mp_cmp(&tb, &ta) != MP_GT) {
			if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
				((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
					goto LBL_ERR;
			}
		}
		if (((res = mp_div_2d(pool, &tb, 1, &tb, NULL)) != MP_OKAY) ||
			((res = mp_div_2d(pool, &tq, 1, &tq, NULL)) != MP_OKAY)) {
			goto LBL_ERR;
		}
	}

/*
	now q == quotient and ta == remainder
 */
	n  = a->sign;
	n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
	if (c != NULL) {
		mp_exch(c, &q);
		c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
	}
	if (d != NULL) {
		mp_exch(d, &ta);
		d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
	}
LBL_ERR:
	_mp_clear_multi(&ta, &tb, &tq, &q, NULL, NULL, NULL, NULL);
	return res;
}
#else /* MP_DIV_SMALL */

int32 mp_div(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c, mp_int * d)
{
	mp_int		q, x, y, t1, t2;
	int32		res, n, t, i, norm, neg;

/*
	is divisor zero ?
 */
	if (mp_iszero(b) == 1) {
		return MP_VAL;
	}

/*
	if a < b then q=0, r = a
 */
	if (mp_cmp_mag(a, b) == MP_LT) {
		if (d != NULL) {
			res = mp_copy(a, d);
		} else {
			res = MP_OKAY;
		}
		if (c != NULL) {
		mp_zero(c);
		}
		return res;
	}

	if ((res = mp_init_size(pool, &q, a->used + 2)) != MP_OKAY) {
		return res;
	}
	q.used = a->used + 2;

	if ((res = mp_init(pool, &t1)) != MP_OKAY) {
		goto LBL_Q;
	}

	if ((res = mp_init(pool, &t2)) != MP_OKAY) {
		goto LBL_T1;
	}

	if ((res = mp_init_copy(pool, &x, a)) != MP_OKAY) {
		goto LBL_T2;
	}

	if ((res = mp_init_copy(pool, &y, b)) != MP_OKAY) {
		goto LBL_X;
	}

/*
	fix the sign
 */
	neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
	x.sign = y.sign = MP_ZPOS;

/*
	normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT]
 */
	norm = mp_count_bits(&y) % DIGIT_BIT;
	if (norm < (int32)(DIGIT_BIT-1)) {
		norm = (DIGIT_BIT-1) - norm;
		if ((res = mp_mul_2d(&x, norm, &x)) != MP_OKAY) {
			goto LBL_Y;
		}
		if ((res = mp_mul_2d(&y, norm, &y)) != MP_OKAY) {
			goto LBL_Y;
		}
	} else {
		norm = 0;
	}

/*
	note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4
 */
	n = x.used - 1;
	t = y.used - 1;

/*
	while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} }
 */
	if ((res = mp_lshd(&y, n - t)) != MP_OKAY) { /* y = y*b**{n-t} */
		goto LBL_Y;
	}

	while (mp_cmp(&x, &y) != MP_LT) {
		++(q.dp[n - t]);
		if ((res = mp_sub(&x, &y, &x)) != MP_OKAY) {
		goto LBL_Y;
		}
	}

/*
	reset y by shifting it back down
 */
	mp_rshd(&y, n - t);

/*
	step 3. for i from n down to (t + 1)
 */
	for (i = n; i >= (t + 1); i--) {
		if (i > x.used) {
			continue;
		}

/*
		step 3.1 if xi == yt then set q{i-t-1} to b-1,
		otherwise set q{i-t-1} to (xi*b + x{i-1})/yt
 */
		if (x.dp[i] == y.dp[t]) {
			q.dp[i - t - 1] = ((((mp_digit)1) << DIGIT_BIT) - 1);
		} else {
			mp_word tmp;
			tmp = ((mp_word) x.dp[i]) << ((mp_word) DIGIT_BIT);
			tmp |= ((mp_word) x.dp[i - 1]);
			tmp /= ((mp_word) y.dp[t]);
			if (tmp > (mp_word) MP_MASK) {
				tmp = MP_MASK;
			}
			q.dp[i - t - 1] = (mp_digit) (tmp & (mp_word) (MP_MASK));
		}

/*
		while (q{i-t-1} * (yt * b + y{t-1})) > 
				xi * b**2 + xi-1 * b + xi-2 

			do q{i-t-1} -= 1; 
 */
		q.dp[i - t - 1] = (q.dp[i - t - 1] + 1) & MP_MASK;
		do {
			q.dp[i - t - 1] = (q.dp[i - t - 1] - 1) & MP_MASK;

/*
			find left hand
 */
			mp_zero (&t1);
			t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
			t1.dp[1] = y.dp[t];
			t1.used = 2;
			if ((res = mp_mul_d (&t1, q.dp[i - t - 1], &t1)) != MP_OKAY) {
				goto LBL_Y;
			}

/*
			find right hand
 */
			t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
			t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
			t2.dp[2] = x.dp[i];
			t2.used = 3;
		} while (mp_cmp_mag(&t1, &t2) == MP_GT);

/*
		step 3.3 x = x - q{i-t-1} * y * b**{i-t-1}
 */
		if ((res = mp_mul_d(&y, q.dp[i - t - 1], &t1)) != MP_OKAY) {
			goto LBL_Y;
		}

		if ((res = mp_lshd(&t1, i - t - 1)) != MP_OKAY) {
			goto LBL_Y;
		}

		if ((res = mp_sub(&x, &t1, &x)) != MP_OKAY) {
			goto LBL_Y;
		}

/*
	if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; }
 */
		if (x.sign == MP_NEG) {
			if ((res = mp_copy(&y, &t1)) != MP_OKAY) {
				goto LBL_Y;
			}
		if ((res = mp_lshd (&t1, i - t - 1)) != MP_OKAY) {
			goto LBL_Y;
		}
		if ((res = mp_add (&x, &t1, &x)) != MP_OKAY) {
			goto LBL_Y;
		}

		q.dp[i - t - 1] = (q.dp[i - t - 1] - 1UL) & MP_MASK;
		}
	}

/*
	now q is the quotient and x is the remainder 
	[which we have to normalize] 
 */

/*
	get sign before writing to c
 */
	x.sign = x.used == 0 ? MP_ZPOS : a->sign;

	if (c != NULL) {
		mp_clamp(&q);
		mp_exch(&q, c);
		c->sign = neg;
	}

	if (d != NULL) {
		mp_div_2d(pool, &x, norm, &x, NULL);
		mp_exch(&x, d);
	}

	res = MP_OKAY;

LBL_Y:mp_clear (&y);
LBL_X:mp_clear (&x);
LBL_T2:mp_clear (&t2);
LBL_T1:mp_clear (&t1);
LBL_Q:mp_clear (&q);
	return res;
}
#endif /* MP_DIV_SMALL */

/******************************************************************************/
/*
	multiplies |a| * |b| and only computes upto digs digits of result 
	HAC pp. 595, Algorithm 14.12  Modified so you can control how many digits
	of output are created.
 */
#ifdef USE_SMALL_WORD
int32 s_mp_mul_digs(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c, int32 digs)
{
	mp_int		t;
	int32			res, pa, pb, ix, iy;
	mp_digit	u;
	mp_word		r;
	mp_digit	tmpx, *tmpt, *tmpy;

/*
	Can we use the fast multiplier?
 */
	if (((digs) < MP_WARRAY) &&
		MIN (a->used, b->used) < 
		(1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
			return fast_s_mp_mul_digs (pool, a, b, c, digs);
		}

		if ((res = mp_init_size(pool, &t, digs)) != MP_OKAY) {
			return res;
		}
		t.used = digs;

/*
		Compute the digits of the product directly
 */
		pa = a->used;
		for (ix = 0; ix < pa; ix++) {
			/* set the carry to zero */
			u = 0;

/*
			Limit ourselves to making digs digits of output.
*/
			pb = MIN (b->used, digs - ix);

/*
			Setup some aliases. Copy of the digit from a used
			within the nested loop
 */
			tmpx = a->dp[ix];

/*
			An alias for the destination shifted ix places
 */
			tmpt = t.dp + ix;

/*
			An alias for the digits of b
 */
			tmpy = b->dp;

/*
			Compute the columns of the output and propagate the carry
 */
			for (iy = 0; iy < pb; iy++) {
				/* compute the column as a mp_word */
				r       = ((mp_word)*tmpt) +
					((mp_word)tmpx) * ((mp_word)*tmpy++) +
					((mp_word) u);

				/* the new column is the lower part of the result */
				*tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));

				/* get the carry word from the result */
				u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
			}
/*
			Set carry if it is placed below digs
 */
			if (ix + iy < digs) {
				*tmpt = u;
			}
		}

		mp_clamp (&t);
		mp_exch (&t, c);

		mp_clear (&t);
		return MP_OKAY;
}
#endif /* USE_SMALL_WORD */

/******************************************************************************/
/*
	Fast (comba) multiplier

	This is the fast column-array [comba] multiplier.  It is designed to
	compute the columns of the product first then handle the carries afterwards.
	This has the effect of making the nested loops that compute the columns
	very simple and schedulable on super-scalar processors.

	This has been modified to produce a variable number of digits of output so
	if say only a half-product is required you don't have to compute the upper
	half (a feature required for fast Barrett reduction).

	Based on Algorithm 14.12 on pp.595 of HAC.
*/

int32 fast_s_mp_mul_digs(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c,
						  int32 digs)
{
	int32		olduse, res, pa, ix, iz, neg;
	mp_digit	W[MP_WARRAY];
	register	mp_word  _W;

	neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;

/*
	grow the destination as required
 */
	if (c->alloc < digs) {
		if ((res = mp_grow(c, digs)) != MP_OKAY) {
			return res;
		}
	}

/*
	number of output digits to produce
 */
	pa = MIN(digs, a->used + b->used);

/*
	clear the carry
 */
	_W = 0;
	for (ix = 0; ix < pa; ix++) { 
		int32		tx, ty;
		int32		iy;
		mp_digit	*tmpx, *tmpy;

/*
		get offsets into the two bignums
 */
		ty = MIN(b->used-1, ix);
		tx = ix - ty;

/*
		setup temp aliases
 */
		tmpx = a->dp + tx;
		tmpy = b->dp + ty;

/*
		this is the number of times the loop will iterrate, essentially its
		while (tx++ < a->used && ty-- >= 0) { ... }
 */
		iy = MIN(a->used-tx, ty+1);

/*
		execute loop
 */
		for (iz = 0; iz < iy; ++iz) {
			_W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
		}

/*
		store term
 */
		W[ix] = (mp_digit)(_W & MP_MASK);

/*
		make next carry
 */
		_W = _W >> ((mp_word)DIGIT_BIT);
	}

/*
	store final carry
 */
	W[ix] = (mp_digit)(_W & MP_MASK);

/*
	setup dest
 */
	olduse  = c->used;
	c->used = pa;

	{
		register mp_digit *tmpc;
		tmpc = c->dp;
		for (ix = 0; ix < pa+1; ix++) {
/*
			now extract the previous digit [below the carry]
 */
			*tmpc++ = W[ix];
		}

/*
		clear unused digits [that existed in the old copy of c]
 */
		for (; ix < olduse; ix++) {
			*tmpc++ = 0;
		}
	}
	mp_clamp (c);
	c->sign = (c->used > 0) ? neg : MP_ZPOS;
	return MP_OKAY;
}

/******************************************************************************/
/*
	b = a*2
 */