mpi.c

PeerSec Networks MatrixSSL?is an embedded SSL implementation designed for small footprint applicatio

	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
 */
int32 mp_mul_2 (mp_int * a, mp_int * b)
{
	int32		x, res, oldused;

/*
	grow to accomodate result
 */
	if (b->alloc < a->used + 1) {
		if ((res = mp_grow (b, a->used + 1)) != MP_OKAY) {
			return res;
		}
	}

	oldused = b->used;
	b->used = a->used;

	{
		register mp_digit r, rr, *tmpa, *tmpb;

		/* alias for source */
		tmpa = a->dp;

		/* alias for dest */
		tmpb = b->dp;

		/* carry */
		r = 0;
		for (x = 0; x < a->used; x++) {

/*
			get what will be the *next* carry bit from the MSB of the 
			current digit 
 */
			rr = *tmpa >> ((mp_digit)(DIGIT_BIT - 1));

/*
			now shift up this digit, add in the carry [from the previous]
 */
			*tmpb++ = ((*tmpa++ << ((mp_digit)1)) | r) & MP_MASK;

/*			copy the carry that would be from the source digit into the next
			iteration 
 */
			r = rr;
		}

/*
		new leading digit?
 */
		if (r != 0) {
/*
			add a MSB which is always 1 at this point
 */
			*tmpb = 1;
			++(b->used);
		}

/*
		now zero any excess digits on the destination that we didn't write to
 */
		tmpb = b->dp + b->used;
		for (x = b->used; x < oldused; x++) {
			*tmpb++ = 0;
		}
	}
	b->sign = a->sign;
	return MP_OKAY;
}

/******************************************************************************/
/*
	multiply by a digit
 */
int32 mp_mul_d(mp_int * a, mp_digit b, mp_int * c)
{
	mp_digit	u, *tmpa, *tmpc;
	mp_word		r;
	int32			ix, res, olduse;

/*
	make sure c is big enough to hold a*b
 */
	if (c->alloc < a->used + 1) {
		if ((res = mp_grow (c, a->used + 1)) != MP_OKAY) {
			return res;
		}
	}

/*
	get the original destinations used count
 */
	olduse = c->used;

/*
	set the sign
 */
	c->sign = a->sign;

/*
	alias for a->dp [source]
 */
	tmpa = a->dp;

/*
	alias for c->dp [dest]
 */
	tmpc = c->dp;

	/* zero carry */
	u = 0;

	/* compute columns */
	for (ix = 0; ix < a->used; ix++) {
/*
		compute product and carry sum for this term
 */
		r       = ((mp_word) u) + ((mp_word)*tmpa++) * ((mp_word)b);

/*
		mask off higher bits to get a single digit
 */
		*tmpc++ = (mp_digit) (r & ((mp_word) MP_MASK));

/*
		send carry into next iteration
 */
		u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
	}

/*
	store final carry [if any] and increment ix offset
 */
	*tmpc++ = u;
	++ix;

/*
	now zero digits above the top
 */
	while (ix++ < olduse) {
		*tmpc++ = 0;
	}

	/* set used count */
	c->used = a->used + 1;
	mp_clamp(c);

	return MP_OKAY;
}

/******************************************************************************/
/*
	low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16
 */
#ifdef USE_SMALL_WORD
int32 s_mp_sqr (psPool_t *pool, mp_int * a, mp_int * b)
{
	mp_int		t;
	int32			res, ix, iy, pa;
	mp_word		r;
	mp_digit	u, tmpx, *tmpt;

	pa = a->used;
	if ((res = mp_init_size(pool, &t, 2*pa + 1)) != MP_OKAY) {
		return res;
	}
	
/*
	default used is maximum possible size
 */
	t.used = 2*pa + 1;

	for (ix = 0; ix < pa; ix++) {