mpi.c

matrix ssl代码

		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[...]
 */
	{
		register mp_word		*_W;
		register mp_digit	*tmpx;

/*		
		Alias for the W[] array
 */
		_W   = W;

/*
		Alias for the digits of  x
 */
		tmpx = x->dp;

/*
		Copy the digits of a into W[0..a->used-1]
 */
		for (ix = 0; ix < x->used; ix++) {
			*_W++ = *tmpx++;
		}

/*
		Zero the high words of W[a->used..m->used*2]
 */
		for (; ix < n->used * 2 + 1; ix++) {
			*_W++ = 0;
		}
	}

/*
	Now we proceed to zero successive digits from the least
	significant upwards.
 */
	for (ix = 0; ix < n->used; ix++) {
/*
		mu = ai * m' mod b

		We avoid a double precision multiplication (which isn't required) by 
		casting the value down to a mp_digit.  Note this requires that
		W[ix-1] have  the carry cleared (see after the inner loop)
 */
		register mp_digit mu;
		mu = (mp_digit) (((W[ix] & MP_MASK) * rho) & MP_MASK);

/*		
		a = a + mu * m * b**i

		This is computed in place and on the fly.  The multiplication by b**i 
		is handled by offseting which columns the results are added to.
		
		Note the comba method normally doesn't handle carries in the inner loop
		In this case we fix the carry from the previous column since the
		Montgomery reduction requires digits of the result (so far) [see above]
		to work.  This is handled by fixing up one carry after the inner loop.
		The carry fixups are done in order so after these loops the first
		m->used words of W[] have the carries fixed
 */
		{
			register int32		iy;
			register mp_digit	*tmpn;
			register mp_word	*_W;

/*
			Alias for the digits of the modulus
 */
			tmpn = n->dp;

/*
			Alias for the columns set by an offset of ix
 */
			_W = W + ix;

/*
				inner loop
 */
				for (iy = 0; iy < n->used; iy++) {
					*_W++ += ((mp_word)mu) * ((mp_word)*tmpn++);
				}
			}

/*
		Now fix carry for next digit, W[ix+1]
 */
		W[ix + 1] += W[ix] >> ((mp_word) DIGIT_BIT);
	}

/*
		Now we have to propagate the carries and shift the words downward [all those 
		least significant digits we zeroed].
 */
	{
		register mp_digit	*tmpx;
		register mp_word	*_W, *_W1;

/*
		Now fix rest of carries 
 */

/*
		alias for current word
 */
		_W1 = W + ix;

/*
		alias for next word, where the carry goes
 */
		_W = W + ++ix;

		for (; ix <= n->used * 2 + 1; ix++) {
			*_W++ += *_W1++ >> ((mp_word) DIGIT_BIT);
		}

/*
		copy out, A = A/b**n

		The result is A/b**n but instead of converting from an
		array of mp_word to mp_digit than calling mp_rshd
		we just copy them in the right order
 */

/*
		alias for destination word
 */
		tmpx = x->dp;

/*
		alias for shifted double precision result
 */
		_W = W + n->used;

		for (ix = 0; ix < n->used + 1; ix++) {
		*tmpx++ = (mp_digit)(*_W++ & ((mp_word) MP_MASK));
		}

/*
		zero oldused digits, if the input a was larger than
		m->used+1 we'll have to clear the digits
 */
		for (; ix < olduse; ix++) {
			*tmpx++ = 0;
		}
	}

/*
	Set the max used and clamp
 */
	x->used = n->used + 1;
	mp_clamp(x);

/*
	if A >= m then A = A - m
 */
	if (mp_cmp_mag(x, n) != MP_LT) {
		return s_mp_sub(x, n, x);
	}
	return MP_OKAY;
}

/******************************************************************************/
/*
	High level addition (handles signs)
 */
int32 mp_add (mp_int * a, mp_int * b, mp_int * c)
{
	int32		sa, sb, res;

/*
	Get sign of both inputs
 */
	sa = a->sign;
	sb = b->sign;

/*
	Handle two cases, not four.
 */
	if (sa == sb) {
/*
		Both positive or both negative. Add their magnitudes, copy the sign.
 */
		c->sign = sa;
		res = s_mp_add (a, b, c);
	} else {
/*
		One positive, the other negative.  Subtract the one with the greater
		magnitude from the one of the lesser magnitude.  The result gets the sign of
		the one with the greater magnitude.
 */
		if (mp_cmp_mag (a, b) == MP_LT) {
			c->sign = sb;
			res = s_mp_sub (b, a, c);
		} else {
			c->sign = sa;
			res = s_mp_sub (a, b, c);
		}
	}
	return res;
}

/******************************************************************************/
/*
	Compare a digit.
 */
int32 mp_cmp_d (mp_int * a, mp_digit b)
{
/*
	Compare based on sign
 */
	if (a->sign == MP_NEG) {
		return MP_LT;
	}

/*
	Compare based on magnitude
 */
	if (a->used > 1) {
		return MP_GT;
	}

/*
	Compare the only digit of a to b
 */
	if (a->dp[0] > b) {
		return MP_GT;
	} else if (a->dp[0] < b) {
		return MP_LT;
	} else {
		return MP_EQ;
	}
}

/******************************************************************************/
/*
	b = a/2
 */
int32 mp_div_2 (mp_int * a, mp_int * b)
{
	int32		x, res, oldused;

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

	oldused = b->used;
	b->used = a->used;
	{
		register mp_digit r, rr, *tmpa, *tmpb;

/*
		Source alias
 */
		tmpa = a->dp + b->used - 1;

/*
		dest alias
 */
		tmpb = b->dp + b->used - 1;

/*
		carry
 */
		r = 0;
		for (x = b->used - 1; x >= 0; x--) {
/*
			Get the carry for the next iteration
 */
			rr = *tmpa & 1;

/*			
			Shift the current digit, add in carry and store
 */
			*tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
/*
			Forward carry to next iteration
 */
			r = rr;
		}

/* 
		Zero excess digits
 */
		tmpb = b->dp + b->used;
		for (x = b->used; x < oldused; x++) {
			*tmpb++ = 0;
		}
	}
	b->sign = a->sign;
	mp_clamp (b);
	return MP_OKAY;
}

/******************************************************************************/
/*
	Computes xR**-1 == x (mod N) via Montgomery Reduction
 */
#ifdef USE_SMALL_WORD
int32 mp_montgomery_reduce (mp_int * x, mp_int * n, mp_digit rho)
{
	int32			ix, res, digs;
	mp_digit	mu;

/*	Can the fast reduction [comba] method be used?

	Note that unlike in mul you're safely allowed *less* than the available
	columns [255 per default] since carries are fixed up in the inner loop.
 */
	digs = n->used * 2 + 1;
	if ((digs < MP_WARRAY) && 
		n->used < 
		(1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
			return fast_mp_montgomery_reduce (x, n, rho);
		}

/*
		Grow the input as required.
 */
		if (x->alloc < digs) {
			if ((res = mp_grow (x, digs)) != MP_OKAY) {
				return res;
			}
		}
		x->used = digs;

		for (ix = 0; ix < n->used; ix++) {
/*
			mu = ai * rho mod b

			The value of rho must be precalculated via mp_montgomery_setup()
			such that it equals -1/n0 mod b this allows the following inner
			loop to reduce the input one digit at a time
 */
			mu = (mp_digit)(((mp_word)x->dp[ix]) * ((mp_word)rho) & MP_MASK);

			/* a = a + mu * m * b**i */
			{
				register int32 iy;
				register mp_digit *tmpn, *tmpx, u;
				register mp_word r;

/*
				alias for digits of the modulus
 */
				tmpn = n->dp;

/*
				alias for the digits of x [the input]
 */
				tmpx = x->dp + ix;

/*
				set the carry to zero
 */
				u = 0;

/*
				Multiply and add in place
 */
				for (iy = 0; iy < n->used; iy++) {
					/* compute product and sum */
					r = ((mp_word)mu) * ((mp_word)*tmpn++) +
						((mp_word) u) + ((mp_word) * tmpx);

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

					/* fix digit */
					*tmpx++ = (mp_digit)(r & ((mp_word) MP_MASK));
				}
				/* At this point the ix'th digit of x should be zero */


/*
				propagate carries upwards as required
 */
				while (u) {
					*tmpx		+= u;
					u			= *tmpx >> DIGIT_BIT;
					*tmpx++ &= MP_MASK;
				}
			}
		}

/*
		At this point the n.used'th least significant digits of x are all zero
		which means we can shift x to the right by n.used digits and the 
		residue is unchanged.
*/
		/* x = x/b**n.used */
		mp_clamp(x);
		mp_rshd (x, n->used);

		/* if x >= n then x = x - n */
		if (mp_cmp_mag (x, n) != MP_LT) {
			return s_mp_sub (x, n, x);
		}

		return MP_OKAY;
}
#endif /* USE_SMALL_WORD */

/******************************************************************************/
/*
	Setups the montgomery reduction stuff.
 */
int32 mp_montgomery_setup (mp_int * n, mp_digit * rho)
{
	mp_digit x, b;

/*
	fast inversion mod 2**k
	
	Based on the fact that
	
	XA = 1 (mod 2**n)	=>  (X(2-XA)) A		= 1 (mod 2**2n)
						=>  2*X*A - X*X*A*A	= 1
						=>  2*(1) - (1)		= 1
*/
	b = n->dp[0];

	if ((b & 1) == 0) {
		return MP_VAL;
	}

	x = (((b + 2) & 4) << 1) + b;		/* here x*a==1 mod 2**4 */
	x = (x * (2 - b * x)) & MP_MASK;	/* here x*a==1 mod 2**8 */
#if !defined(MP_8BIT)
	x = (x * (2 - b * x)) & MP_MASK;	/* here x*a==1 mod 2**8 */
#endif /* MP_8BIT */
#if defined(MP_64BIT) || !(defined(MP_8BIT) || defined(MP_16BIT))
	x *= 2 - b * x;						/* here x*a==1 mod 2**32 */
#endif
#ifdef MP_64BIT
	x *= 2 - b * x;						/* here x*a==1 mod 2**64 */
#endif /* MP_64BIT */

	/* rho = -1/m mod b */
	*rho = (((mp_word) 1 << ((mp_word) DIGIT_BIT)) - x) & MP_MASK;

	return MP_OKAY;
}

/******************************************************************************/
/*
	High level subtraction (handles signs)
 */
int32 mp_sub (mp_int * a, mp_int * b, mp_int * c)
{
	int32		sa, sb, res;

	sa = a->sign;
	sb = b->sign;

	if (sa != sb) {
/*
		Subtract a negative from a positive, OR subtract a positive from a
		negative.  In either case, ADD their magnitudes, and use the sign of
		the first number.
 */
		c->sign = sa;
		res = s_mp_add (a, b, c);
	} else {
/*
		Subtract a positive from a positive, OR subtract a negative 
		from a negative. First, take the difference between their
		magnitudes, then...
 */
		if (mp_cmp_mag (a, b) != MP_LT) {
/*
			Copy the sign from the first
 */
			c->sign = sa;
			/* The first has a larger or equal magnitude */
			res = s_mp_sub (a, b, c);
		} else {
/*
			The result has the *opposite* sign from the first number.
 */
			c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
/*
			The second has a larger magnitude 
 */
			res = s_mp_sub (b, a, c);
		}
	}
	return res;
}

/******************************************************************************/
/*
	calc a value mod 2**b
 */
int32 mp_mod_2d (mp_int * a, int32 b, mp_int * c)
{
	int32		x, res;

/*
	if b is <= 0 then zero the int32
 */
	if (b <= 0) {
		mp_zero (c);
		return MP_OKAY;
	}

/*
	If the modulus is larger than the value than return
 */
	if (b >=(int32) (a->used * DIGIT_BIT)) {
		res = mp_copy (a, c);
		return res;
	}

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

/*
	Zero digits above the last digit of the modulus
 */
	for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
		c->dp[x] = 0;
	}
/*
	Clear the digit that is not completely outside/inside the modulus
 */
	c->dp[b / DIGIT_BIT] &=
		(mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1));
	mp_clamp (c);
	return MP_OKAY;
}

/******************************************************************************/
/*
	Shift right a certain amount of digits.
 */
void mp_rshd (mp_int * a, int32 b)
{
	int32		x;

/*
	If b <= 0 then ignore it
 */
	if (b <= 0) {
		return;
	}

/*
	If b > used then simply zero it and return.
*/
	if (a->used <= b) {
		mp_zero (a);