mpi.c

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

/*
		first calculate the digit at 2*ix
		calculate double precision result
 */
		r = ((mp_word) t.dp[2*ix]) +
			((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);

/*
		store lower part in result
 */
		t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));

/*
		get the carry
 */
		u = (mp_digit)(r >> ((mp_word) DIGIT_BIT));

/*
		left hand side of A[ix] * A[iy]
 */
		tmpx = a->dp[ix];

/*
		alias for where to store the results
 */
		tmpt = t.dp + (2*ix + 1);

		for (iy = ix + 1; iy < pa; iy++) {
/*
			first calculate the product
 */
			r = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);

/*
			now calculate the double precision result, note we use addition
			instead of *2 since it's easier to optimize
 */
			r       = ((mp_word) *tmpt) + r + r + ((mp_word) u);

/*
			store lower part
 */
			*tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));

			/* get carry */
			u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
		}
		/* propagate upwards */
		while (u != ((mp_digit) 0)) {
			r       = ((mp_word) *tmpt) + ((mp_word) u);
			*tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
			u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
		}
	}

	mp_clamp (&t);
	mp_exch (&t, b);
	mp_clear (&t);
	return MP_OKAY;
}
#endif /* USE_SMALL_WORD */

/******************************************************************************/
/*
	fast squaring

	This is the comba method where the columns of the product are computed
	first then the carries are computed.  This has the effect of making a very
	simple inner loop that is executed the most

	W2 represents the outer products and W the inner.

	A further optimizations is made because the inner products are of the 
	form "A * B * 2".  The *2 part does not need to be computed until the end
	which is good because 64-bit shifts are slow!

	Based on Algorithm 14.16 on pp.597 of HAC.

	This is the 1.0 version, but no SSE stuff
*/
int32 fast_s_mp_sqr(psPool_t *pool, mp_int * a, mp_int * b)
{
	int32		olduse, res, pa, ix, iz;
	mp_digit	W[MP_WARRAY], *tmpx;
	mp_word		W1;

/* 
	grow the destination as required
 */
	pa = a->used + a->used;
	if (b->alloc < pa) {
		if ((res = mp_grow(b, pa)) != MP_OKAY) {
			return res;
		}
	}

/*
	number of output digits to produce
 */
	W1 = 0;
	for (ix = 0; ix < pa; ix++) { 
		int32		tx, ty, iy;
		mp_word		_W;
		mp_digit	*tmpy;

/*
		clear counter
 */
		_W = 0;

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

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

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

/*
		now for squaring tx can never equal ty 
		we halve the distance since they approach at a rate of 2x
		and we have to round because odd cases need to be executed
*/
		iy = MIN(iy, (ty-tx+1)>>1);

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

/*
		double the inner product and add carry
 */
		_W = _W + _W + W1;

/*
		even columns have the square term in them
 */
		if ((ix&1) == 0) {
			_W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
		}

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

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

/*
	setup dest
 */
	olduse  = b->used;
	b->used = a->used+a->used;

	{
		mp_digit *tmpb;
		tmpb = b->dp;
		for (ix = 0; ix < pa; ix++) {
			*tmpb++ = W[ix] & MP_MASK;
		}

/*
		clear unused digits [that existed in the old copy of c]
 */
		for (; ix < olduse; ix++) {
			*tmpb++ = 0;
		}
	}
	mp_clamp(b);
	return MP_OKAY;
}

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

	Simple algorithm which zeroes the int32, grows it then just sets one bit
	as required.
 */
int32 mp_2expt (mp_int * a, int32 b)
{
	int32		res;

/*
	zero a as per default
 */
	mp_zero (a);

/*
	grow a to accomodate the single bit
 */
	if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
		return res;
	}

/*
	set the used count of where the bit will go
 */
	a->used = b / DIGIT_BIT + 1;

/*
	put the single bit in its place
 */
	a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);

	return MP_OKAY;
}

/******************************************************************************/
/*
	init an mp_init for a given size
 */
int32 mp_init_size(psPool_t *pool, mp_int * a, int32 size)
{
	int		x;
/*
	pad size so there are always extra digits
 */
	size += (MP_PREC * 2) - (size % MP_PREC);	

/*
	alloc mem
 */
	a->dp = OPT_CAST(mp_digit) psMalloc(pool, sizeof (mp_digit) * size);
	if (a->dp == NULL) {
		return MP_MEM;
	}
	a->used  = 0;
	a->alloc = size;
	a->sign  = MP_ZPOS;

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

/******************************************************************************/
/*
	low level addition, based on HAC pp.594, Algorithm 14.7
 */
int32 s_mp_add (mp_int * a, mp_int * b, mp_int * c)
{
	mp_int		*x;
	int32			olduse, res, min, max;

/*
	find sizes, we let |a| <= |b| which means we have to sort them.  "x" will
	point to the input with the most digits
 */
	if (a->used > b->used) {
		min = b->used;
		max = a->used;
		x = a;
	} else {
		min = a->used;
		max = b->used;
		x = b;
	}

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

/*
	get old used digit count and set new one
 */
	olduse = c->used;
	c->used = max + 1;

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

		/* alias for digit pointers */

		/* first input */
		tmpa = a->dp;

		/* second input */
		tmpb = b->dp;

		/* destination */
		tmpc = c->dp;

		/* zero the carry */
		u = 0;
		for (i = 0; i < min; i++) {
/*
			Compute the sum at one digit, T[i] = A[i] + B[i] + U
 */
			*tmpc = *tmpa++ + *tmpb++ + u;

/*
			U = carry bit of T[i]
 */
			u = *tmpc >> ((mp_digit)DIGIT_BIT);

/*
			take away carry bit from T[i]
 */
			*tmpc++ &= MP_MASK;
		}

/*
		now copy higher words if any, that is in A+B if A or B has more digits add
		those in 
 */
		if (min != max) {
			for (; i < max; i++) {
				/* T[i] = X[i] + U */
				*tmpc = x->dp[i] + u;

				/* U = carry bit of T[i] */
				u = *tmpc >> ((mp_digit)DIGIT_BIT);

				/* take away carry bit from T[i] */
				*tmpc++ &= MP_MASK;
			}
		}

		/* add carry */
		*tmpc++ = u;

/*
		clear digits above oldused
 */
		for (i = c->used; i < olduse; i++) {
			*tmpc++ = 0;
		}
	}

	mp_clamp (c);
	return MP_OKAY;
}

/******************************************************************************/
#ifdef USE_SMALL_WORD
/*
	FUTURE - this is never needed, SLOW or not, because RSA exponents are
	always odd. 
*/
int32 mp_invmod(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c)
{
	mp_int		x, y, u, v, A, B, C, D;
	int32			res;

/*
	b cannot be negative
 */
	if (b->sign == MP_NEG || mp_iszero(b) == 1) {
		return MP_VAL;
	}

/*
	if the modulus is odd we can use a faster routine instead
 */
	if (mp_isodd (b) == 1) {
		return fast_mp_invmod(pool, a, b, c);
	}

/*
	init temps
 */
	if ((res = _mp_init_multi(pool, &x, &y, &u, &v,
			&A, &B, &C, &D)) != MP_OKAY) {
		return res;
	}

	/* x = a, y = b */
	if ((res = mp_copy(a, &x)) != MP_OKAY) {
		goto LBL_ERR;
	}
	if ((res = mp_copy(b, &y)) != MP_OKAY) {
		goto LBL_ERR;
	}

/*
	2. [modified] if x,y are both even then return an error!
 */
	if (mp_iseven(&x) == 1 && mp_iseven (&y) == 1) {
		res = MP_VAL;
		goto LBL_ERR;
	}

/*
	3. u=x, v=y, A=1, B=0, C=0,D=1
 */
	if ((res = mp_copy(&x, &u)) != MP_OKAY) {
		goto LBL_ERR;
	}
	if ((res = mp_copy(&y, &v)) != MP_OKAY) {
		goto LBL_ERR;
	}
	mp_set (&A, 1);
	mp_set (&D, 1);

top:
/*
	4.  while u is even do
 */
	while (mp_iseven(&u) == 1) {
		/* 4.1 u = u/2 */
		if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
		goto LBL_ERR;
		}
		/* 4.2 if A or B is odd then */
		if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
			/* A = (A+y)/2, B = (B-x)/2 */
			if ((res = mp_add(&A, &y, &A)) != MP_OKAY) {
				goto LBL_ERR;
			}
			if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
				goto LBL_ERR;
			}
		}
		/* A = A/2, B = B/2 */
		if ((res = mp_div_2(&A, &A)) != MP_OKAY) {
			goto LBL_ERR;
		}
		if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	5.  while v is even do
 */
	while (mp_iseven(&v) == 1) {
		/* 5.1 v = v/2 */
		if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
		goto LBL_ERR;
		}
		/* 5.2 if C or D is odd then */
		if (mp_isodd(&C) == 1 || mp_isodd (&D) == 1) {
			/* C = (C+y)/2, D = (D-x)/2 */
			if ((res = mp_add(&C, &y, &C)) != MP_OKAY) {
				goto LBL_ERR;
			}
			if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
				goto LBL_ERR;
			}
		}
		/* C = C/2, D = D/2 */
		if ((res = mp_div_2(&C, &C)) != MP_OKAY) {
			goto LBL_ERR;
		}
		if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	6.  if u >= v then
 */
	if (mp_cmp(&u, &v) != MP_LT) {
		/* u = u - v, A = A - C, B = B - D */
		if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
		goto LBL_ERR;
		}

		if ((res = mp_sub(&A, &C, &A)) != MP_OKAY) {
		goto LBL_ERR;
		}

		if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
		goto LBL_ERR;
		}
	} else {
		/* v - v - u, C = C - A, D = D - B */
		if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
		goto LBL_ERR;
		}

		if ((res = mp_sub(&C, &A, &C)) != MP_OKAY) {
		goto LBL_ERR;
		}

		if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
		goto LBL_ERR;
		}
	}

/*
	if not zero goto step 4
 */
	if (mp_iszero(&u) == 0)
		goto top;

/*
	now a = C, b = D, gcd == g*v
 */

/*
	if v != 1 then there is no inverse
 */
	if (mp_cmp_d(&v, 1) != MP_EQ) {
		res = MP_VAL;
		goto LBL_ERR;
	}

/*
	if its too low
 */
	while (mp_cmp_d(&C, 0) == MP_LT) {
		if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	too big
 */
	while (mp_cmp_mag(&C, b) != MP_LT) {
		if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	C is now the inverse
 */
	mp_exch(&C, c);
	res = MP_OKAY;
LBL_ERR:_mp_clear_multi(&x, &y, &u, &v, &A, &B, &C, &D);
	return res;
}
#endif /* USE_SMALL_WORD */

/******************************************************************************/

/*
 *	Computes the modular inverse via binary extended euclidean algorithm,
 *	that is c = 1/a mod b 
 *
 * Based on slow invmod except this is optimized for the case where b is 
 * odd as per HAC Note 14.64 on pp. 610
 */
int32 fast_mp_invmod(psPool_t *pool, mp_int * a, mp_int * b, mp_int * c)
{
	mp_int		x, y, u, v, B, D;
	int32			res, neg;

/*
	2. [modified] b must be odd
 */
	if (mp_iseven (b) == 1) {
		return MP_VAL;
	}

/*
	init all our temps
 */
	if ((res = _mp_init_multi(pool, &x, &y, &u, &v, &B, &D, NULL, NULL)) != MP_OKAY) {
		return res;
	}

/*
	x == modulus, y == value to invert
 */
	if ((res = mp_copy(b, &x)) != MP_OKAY) {
		goto LBL_ERR;
	}

/*
	we need y = |a|
 */
	if ((res = mp_mod(pool, a, b, &y)) != MP_OKAY) {
		goto LBL_ERR;
	}

/*
	3. u=x, v=y, A=1, B=0, C=0,D=1
 */
	if ((res = mp_copy(&x, &u)) != MP_OKAY) {
		goto LBL_ERR;
	}
	if ((res = mp_copy(&y, &v)) != MP_OKAY) {
		goto LBL_ERR;
	}
	mp_set(&D, 1);

top:
/*
	4.  while u is even do
*/
	while (mp_iseven(&u) == 1) {
		/* 4.1 u = u/2 */
		if ((res = mp_div_2(&u, &u)) != MP_OKAY) {
			goto LBL_ERR;
		}
		/* 4.2 if B is odd then */
		if (mp_isodd(&B) == 1) {
			if ((res = mp_sub(&B, &x, &B)) != MP_OKAY) {
				goto LBL_ERR;
			}
		}
		/* B = B/2 */
		if ((res = mp_div_2(&B, &B)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	5.  while v is even do
 */
	while (mp_iseven(&v) == 1) {
		/* 5.1 v = v/2 */
		if ((res = mp_div_2(&v, &v)) != MP_OKAY) {
			goto LBL_ERR;
		}
		/* 5.2 if D is odd then */
		if (mp_isodd(&D) == 1) {
			/* D = (D-x)/2 */
			if ((res = mp_sub(&D, &x, &D)) != MP_OKAY) {
				goto LBL_ERR;
			}
		}
		/* D = D/2 */
		if ((res = mp_div_2(&D, &D)) != MP_OKAY) {
			goto LBL_ERR;
		}
	}

/*
	6.  if u >= v then
 */
	if (mp_cmp(&u, &v) != MP_LT) {
		/* u = u - v, B = B - D */
		if ((res = mp_sub(&u, &v, &u)) != MP_OKAY) {
			goto LBL_ERR;
		}

		if ((res = mp_sub(&B, &D, &B)) != MP_OKAY) {
			goto LBL_ERR;
		}
	} else {
		/* v - v - u, D = D - B */
		if ((res = mp_sub(&v, &u, &v)) != MP_OKAY) {
			goto LBL_ERR;
		}

		if ((res = mp_sub(&D, &B, &D)) != MP_OKAY) {
		goto LBL_ERR;
		}
	}

/*
	if not zero goto step 4
 */
	if (mp_iszero(&u) == 0) {
		goto top;
	}

/*
	now a = C, b = D, gcd == g*v
 */

/*
	if v != 1 then there is no inverse
 */
	if (mp_cmp_d(&v, 1) !=