mpi.c

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

	buf    = 0;
	digidx = X->used - 1;
	bitcpy = 0;
	bitbuf = 0;

	for (;;) {
/*
		grab next digit as required
 */
		if (--bitcnt == 0) {
			/* if digidx == -1 we are out of digits so break */
			if (digidx == -1) {
				break;
			}
			/* read next digit and reset bitcnt */
			buf    = X->dp[digidx--];
			bitcnt = (int)DIGIT_BIT;
		}

		/* grab the next msb from the exponent */
		y     = (mp_digit)(buf >> (DIGIT_BIT - 1)) & 1;
		buf <<= (mp_digit)1;

/*
		if the bit is zero and mode == 0 then we ignore it
		These represent the leading zero bits before the first 1 bit
		in the exponent.  Technically this opt is not required but it
		does lower the # of trivial squaring/reductions used
*/
		if (mode == 0 && y == 0) {
			continue;
		}

/*
		if the bit is zero and mode == 1 then we square
 */
		if (mode == 1 && y == 0) {
			if ((err = mp_sqr (pool, &res, &res)) != MP_OKAY) {
				goto LBL_RES;
			}
			if ((err = redux (&res, P, mp)) != MP_OKAY) {
				goto LBL_RES;
			}
			continue;
		}

/*
		else we add it to the window
 */
		bitbuf |= (y << (winsize - ++bitcpy));
		mode    = 2;

		if (bitcpy == winsize) {
/*
			ok window is filled so square as required and multiply
			square first
 */
			for (x = 0; x < winsize; x++) {
				if ((err = mp_sqr(pool, &res, &res)) != MP_OKAY) {
					goto LBL_RES;
				}
				if ((err = redux(&res, P, mp)) != MP_OKAY) {
					goto LBL_RES;
				}
			}

			/* then multiply */
			if ((err = mp_mul(pool, &res, &M[bitbuf], &res)) != MP_OKAY) {
				goto LBL_RES;
			}
			if ((err = redux(&res, P, mp)) != MP_OKAY) {
				goto LBL_RES;
			}

/*
			empty window and reset
 */
			bitcpy = 0;
			bitbuf = 0;
			mode   = 1;
		}
	}

/*
	if bits remain then square/multiply
 */
	if (mode == 2 && bitcpy > 0) {
		/* square then multiply if the bit is set */
		for (x = 0; x < bitcpy; x++) {
			if ((err = mp_sqr(pool, &res, &res)) != MP_OKAY) {
				goto LBL_RES;
			}
			if ((err = redux(&res, P, mp)) != MP_OKAY) {
				goto LBL_RES;
			}

/*
			get next bit of the window
 */
			bitbuf <<= 1;
			if ((bitbuf & (1 << winsize)) != 0) {
/*
				then multiply
 */
				if ((err = mp_mul(pool, &res, &M[1], &res)) != MP_OKAY) {
					goto LBL_RES;
				}
				if ((err = redux(&res, P, mp)) != MP_OKAY) {
					goto LBL_RES;
				}
			}
		}
	}

/*
	fixup result if Montgomery reduction is used
	recall that any value in a Montgomery system is
	actually multiplied by R mod n.  So we have
	to reduce one more time to cancel out the factor of R.
*/
	if ((err = redux(&res, P, mp)) != MP_OKAY) {
		goto LBL_RES;
	}

/*
	swap res with Y
 */
	mp_exch(&res, Y);
	err = MP_OKAY;
LBL_RES:mp_clear(&res);
LBL_M:
	mp_clear(&M[1]);
	for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
		mp_clear(&M[x]);
	}
	return err;
}

/******************************************************************************/
/*
	Grow as required
 */
int32 mp_grow (mp_int * a, int32 size)
{
	int32			i;
	mp_digit	*tmp;

/* 
	If the alloc size is smaller alloc more ram.
 */
	if (a->alloc < size) {
/*
		ensure there are always at least MP_PREC digits extra on top
 */
		size += (MP_PREC * 2) - (size % MP_PREC);

/*
		Reallocate the array a->dp

		We store the return in a temporary variable in case the operation 
		failed we don't want to overwrite the dp member of a.
*/
		tmp = OPT_CAST(mp_digit) psRealloc(a->dp, sizeof (mp_digit) * size);
		if (tmp == NULL) {
/*
			reallocation failed but "a" is still valid [can be freed]
 */
			return MP_MEM;
		}

/*
		reallocation succeeded so set a->dp
 */
		a->dp = tmp;

/*
		zero excess digits
 */
		i			= a->alloc;
		a->alloc	= size;
		for (; i < a->alloc; i++) {
			a->dp[i] = 0;
		}
	}
	return MP_OKAY;
}

/******************************************************************************/
/*
	b = |a|

	Simple function copies the input and fixes the sign to positive
*/
int32 mp_abs (mp_int * a, mp_int * b)
{
	int32		res;

/*
	copy a to b
 */
	if (a != b) {
		if ((res = mp_copy (a, b)) != MP_OKAY) {
			return res;
		}
	}

/*
	Force the sign of b to positive
 */
	b->sign = MP_ZPOS;

	return MP_OKAY;
}

/******************************************************************************/
/*
	Creates "a" then copies b into it
 */
int32 mp_init_copy(psPool_t *pool, mp_int * a, mp_int * b)
{
	int32		res;

	if ((res = mp_init(pool, a)) != MP_OKAY) {
		return res;
	}
	return mp_copy (b, a);
}

/******************************************************************************/
/* 
	Reverse an array, used for radix code
 */
void bn_reverse (unsigned char *s, int32 len)
{
	int32				ix, iy;
	unsigned char	t;

	ix = 0;
	iy = len - 1;
	while (ix < iy) {
		t		= s[ix];
		s[ix]	= s[iy];
		s[iy]	= t;
		++ix;
		--iy;
	}
}

/******************************************************************************/
/*
	Shift right by a certain bit count (store quotient in c, optional 
	remainder in d)
 */
int32 mp_div_2d(psPool_t *pool, mp_int * a, int32 b, mp_int * c, mp_int * d)
{
	mp_digit	D, r, rr;
	int32			x, res;
	mp_int		t;

/*
	If the shift count is <= 0 then we do no work
 */
	if (b <= 0) {
		res = mp_copy (a, c);
		if (d != NULL) {
			mp_zero (d);
		}
		return res;
	}

	if ((res = mp_init(pool, &t)) != MP_OKAY) {
		return res;
	}

/*
	Get the remainder
 */
	if (d != NULL) {
		if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
			mp_clear (&t);
			return res;
		}
	}

	/* copy */
	if ((res = mp_copy (a, c)) != MP_OKAY) {
		mp_clear (&t);
		return res;
	}

/*
	Shift by as many digits in the bit count
 */
	if (b >= (int32)DIGIT_BIT) {
		mp_rshd (c, b / DIGIT_BIT);
	}

	/* shift any bit count < DIGIT_BIT */
	D = (mp_digit) (b % DIGIT_BIT);
	if (D != 0) {
		register mp_digit *tmpc, mask, shift;

		/* mask */
		mask = (((mp_digit)1) << D) - 1;

		/* shift for lsb */
		shift = DIGIT_BIT - D;

		/* alias */
		tmpc = c->dp + (c->used - 1);

		/* carry */
		r = 0;
		for (x = c->used - 1; x >= 0; x--) {
/*	
			Get the lower  bits of this word in a temp.
 */
			rr = *tmpc & mask;

/*
			shift the current word and mix in the carry bits from the previous word
 */
			*tmpc = (*tmpc >> D) | (r << shift);
			--tmpc;

/*
			set the carry to the carry bits of the current word found above
 */
			r = rr;
		}
	}
	mp_clamp (c);
	if (d != NULL) {
		mp_exch (&t, d);
	}
	mp_clear (&t);
	return MP_OKAY;
}

/******************************************************************************/
/*
	copy, b = a
 */
int32 mp_copy (mp_int * a, mp_int * b)
{
	int32		res, n;

/*
	If dst == src do nothing
 */
	if (a == b) {
		return MP_OKAY;
	}

/*
	Grow dest
 */
	if (b->alloc < a->used) {
		if ((res = mp_grow (b, a->used)) != MP_OKAY) {
			return res;
		}
	}

/*
	Zero b and copy the parameters over
 */
	{
		register mp_digit *tmpa, *tmpb;

		/* pointer aliases */
		/* source */
		tmpa = a->dp;

		/* destination */
		tmpb = b->dp;

		/* copy all the digits */
		for (n = 0; n < a->used; n++) {
			*tmpb++ = *tmpa++;
		}

		/* clear high digits */
		for (; n < b->used; n++) {
			*tmpb++ = 0;
		}
	}

/*
	copy used count and sign
 */
	b->used = a->used;
	b->sign = a->sign;
	return MP_OKAY;
}

/******************************************************************************/
/*
	Returns the number of bits in an int32
 */
int32 mp_count_bits (mp_int * a)
{
	int32			r;
	mp_digit	q;

/* 
	Shortcut
 */
	if (a->used == 0) {
		return 0;
	}

/*
	Get number of digits and add that.
 */
	r = (a->used - 1) * DIGIT_BIT;

/*
	Take the last digit and count the bits in it.
 */
	q = a->dp[a->used - 1];
	while (q > ((mp_digit) 0)) {
		++r;
		q >>= ((mp_digit) 1);
	}
	return r;
}

/******************************************************************************/
/*
	Shift left a certain amount of digits.
 */
int32 mp_lshd (mp_int * a, int32 b)
{
	int32		x, res;

/*
	If its less than zero return.
 */
	if (b <= 0) {
		return MP_OKAY;
	}

/*
	Grow to fit the new digits.
 */
	if (a->alloc < a->used + b) {
		if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
			return res;
		}
	}

	{
		register mp_digit *top, *bottom;

/*
		Increment the used by the shift amount then copy upwards.
 */
		a->used += b;

		/* top */
		top = a->dp + a->used - 1;

		/* base */
		bottom = a->dp + a->used - 1 - b;

/*
		Much like mp_rshd this is implemented using a sliding window
		except the window goes the otherway around.  Copying from
		the bottom to the top.  see bn_mp_rshd.c for more info.
 */
		for (x = a->used - 1; x >= b; x--) {
			*top-- = *bottom--;
		}

		/* zero the lower digits */
		top = a->dp;
		for (x = 0; x < b; x++) {
			*top++ = 0;
		}
	}
	return MP_OKAY;
}

/******************************************************************************/
/*
	Set to a digit.
 */
void mp_set (mp_int * a, mp_digit b)
{
	mp_zero (a);
	a->dp[0] = b & MP_MASK;
	a->used  = (a->dp[0] != 0) ? 1 : 0;
}

/******************************************************************************/
/*
	Swap the elements of two integers, for cases where you can't simply swap 
	the 	mp_int pointers around 
*/
void mp_exch (mp_int * a, mp_int * b)
{
	mp_int		t;

	t	= *a;
	*a	= *b;
	*b	= t;
}

/******************************************************************************/
/*
	High level multiplication (handles sign)
 */
int32 mp_mul(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c)
{
	int32			res, neg;
	int32			digs = a->used + b->used + 1;
	
	neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;

/*	Can we use the fast multiplier?
	
	The fast multiplier can be used if the output will have less than 
	MP_WARRAY digits and the number of digits won't affect carry propagation
*/
	if ((digs < MP_WARRAY) && MIN(a->used, b->used) <= 
			(1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
		res = fast_s_mp_mul_digs(pool, a, b, c, digs);
	} else {
		res = s_mp_mul(pool, a, b, c);
	}
	c->sign = (c->used > 0) ? neg : MP_ZPOS;
	return res;
}

/******************************************************************************/
/* 
	c = a mod b, 0 <= c < b
 */
int32 mp_mod(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c)
{
	mp_int		t;
	int32			res;

	if ((res = mp_init(pool, &t)) != MP_OKAY) {
		return res;
	}

	if ((res = mp_div (pool, a, b, NULL, &t)) != MP_OKAY) {
		mp_clear (&t);
		return res;
	}

	if (t.sign != b->sign) {
		res = mp_add (b, &t, c);
	} else {
		res = MP_OKAY;
		mp_exch (&t, c);
	}

	mp_clear (&t);
	return res;
}

/******************************************************************************/
/*
	shifts with subtractions when the result is greater than b.

	The method is slightly modified to shift B unconditionally upto just under
	the leading bit of b.  This saves alot of multiple precision shifting.
*/
int32 mp_montgomery_calc_normalization (mp_int * a, mp_int * b)
{
	int32		x, bits, res;

/*
	How many bits of last digit does b use
 */
	bits = mp_count_bits (b) % DIGIT_BIT;

	if (b->used > 1) {
		if ((res = mp_2expt(a, (b->used - 1) * DIGIT_BIT + bits - 1)) != MP_OKAY) {
			return res;
		}
	} else {
		mp_set(a, 1);
		bits = 1;
	}

/*
	Now compute C = A * B mod b
 */
	for (x = bits - 1; x < (int32)DIGIT_BIT; x++) {
		if ((res = mp_mul_2(a, a)) != MP_OKAY) {
			return res;
		}
		if (mp_cmp_mag(a, b) != MP_LT) {
			if ((res = s_mp_sub(a, b, a)) != MP_OKAY) {
				return res;
			}
		}
	}

	return MP_OKAY;
}

/******************************************************************************/
/*
	d = a * b (mod c)
 */
int32 mp_mulmod(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c, mp_int * d)
{
	int32		res;
	mp_int		t;

	if ((res = mp_init(pool, &t)) != MP_OKAY) {
		return res;
	}

	if ((res = mp_mul (pool, a, b, &t)) != MP_OKAY) {
		mp_clear (&t);
		return res;
	}
	res = mp_mod (pool, &t, c, d);
	mp_clear (&t);
	return res;
}

/******************************************************************************/
/*
	Computes b = a*a
 */
#ifdef USE_SMALL_WORD
int32 mp_sqr (psPool_t *pool, mp_int * a, mp_int * b)
{
	int32		res;

/*
	Can we use the fast comba multiplier?
 */
	if ((a->used * 2 + 1) < MP_WARRAY && a->used < 
			(1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
		res = fast_s_mp_sqr (pool, a, b);
	} else {
		res = s_mp_sqr (pool, a, b);
	}
	b->sign = MP_ZPOS;
	return res;
}
#endif /* USE_SMALL_WORD */

/******************************************************************************/
/* 
	Computes xR**-1 == x (mod N) via Montgomery Reduction.

	This is an optimized implementation of montgomery_reduce 
	which uses the comba method to quickly calculate the columns of the
	reduction.

	Based on Algorithm 14.32 on pp.601 of HAC.
*/

int32 fast_mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
{
	int32		ix, res, olduse;
/*
	FUTURE - lower this stack usage, this is around 1K!.
 */
	mp_word		W[MP_WARRAY];

/*
	Get old used count
 */
	olduse = x->used;

/*
	Grow a as required
 */
	if (x->alloc < n->used + 1) {
		if ((res = mp_grow(x, n->used + 1)) != MP_OKAY) {
			return res;
		}
	}

/*
	First we have to get the digits of the input into
	an array of double precision words W[...]
 */
	{