patch to Vanilla Tomato 1.28

1 files changed, 3667 insertions, 0 deletions

diff --git a/release/src/router/matrixssl/src/crypto/peersec/mpi.c b/release/src/router/matrixssl/src/crypto/peersec/mpi.c
new file mode 100644
index 00000000..c37353d3
--- /dev/null
+++ b/release/src/router/matrixssl/src/crypto/peersec/mpi.c
@@ -0,0 +1,3667 @@
+/*	
+ *	mpi.c
+ *	Release $Name: MATRIXSSL_1_8_8_OPEN $
+ *
+ *	multiple-precision integer library
+ */
+/*
+ *	Copyright (c) PeerSec Networks, 2002-2009. All Rights Reserved.
+ *	The latest version of this code is available at http://www.matrixssl.org
+ *
+ *	This software is open source; you can redistribute it and/or modify
+ *	it under the terms of the GNU General Public License as published by
+ *	the Free Software Foundation; either version 2 of the License, or
+ *	(at your option) any later version.
+ *
+ *	This General Public License does NOT permit incorporating this software 
+ *	into proprietary programs.  If you are unable to comply with the GPL, a 
+ *	commercial license for this software may be purchased from PeerSec Networks
+ *	at http://www.peersec.com
+ *	
+ *	This program is distributed in WITHOUT ANY WARRANTY; without even the 
+ *	implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
+ *	See the GNU General Public License for more details.
+ *	
+ *	You should have received a copy of the GNU General Public License
+ *	along with this program; if not, write to the Free Software
+ *	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *	http://www.gnu.org/copyleft/gpl.html
+ */
+/******************************************************************************/
+
+#include "../cryptoLayer.h"
+#include <stdarg.h>
+
+#ifndef USE_MPI2
+
+static int32 mp_exptmod_fast (psPool_t *pool, mp_int * G, mp_int * X,
+				mp_int * P, mp_int * Y, int32 redmode);
+
+/******************************************************************************/
+/*
+	FUTURE
+	1. Convert the mp_init and mp_clear functions to not use malloc + free,
+	but to use static storage within the bignum variable instead - but
+	how to handle grow()?  Maybe use a simple memory allocator
+	2. verify stack usage of all functions and use of MP_LOW_MEM:
+		fast_mp_montgomery_reduce
+		fast_s_mp_mul_digs
+		fast_s_mp_sqr
+		fast_s_mp_mul_high_digs
+	3. HAC stands for Handbook of Applied Cryptography
+		http://www.cacr.math.uwaterloo.ca/hac/
+*/
+/******************************************************************************/
+/*
+	Utility functions
+*/
+void psZeromem(void *dst, size_t len)
+{
+	unsigned char *mem = (unsigned char *)dst;
+	
+	if (dst == NULL) {
+		return;
+	}
+	while (len-- > 0) {
+		*mem++ = 0;
+	}
+}
+
+void psBurnStack(unsigned long len)
+{
+	unsigned char buf[32];
+	
+	psZeromem(buf, sizeof(buf));
+	if (len > (unsigned long)sizeof(buf)) {
+		psBurnStack(len - sizeof(buf));
+	}
+}
+
+/******************************************************************************/
+/*
+	Multiple precision integer functions
+	Note: we don't use va_args here to prevent portability issues.
+*/
+int32 _mp_init_multi(psPool_t *pool, mp_int *mp0, mp_int *mp1, mp_int *mp2, 
+					mp_int *mp3, mp_int *mp4, mp_int *mp5, 
+					mp_int *mp6, mp_int *mp7)
+{
+	mp_err	res		= MP_OKAY;		/* Assume ok until proven otherwise */
+	int32		n		= 0;			/* Number of ok inits */
+	mp_int	*tempArray[9];
+	
+	tempArray[0] = mp0;
+	tempArray[1] = mp1;
+	tempArray[2] = mp2;
+	tempArray[3] = mp3;
+	tempArray[4] = mp4;
+	tempArray[5] = mp5;
+	tempArray[6] = mp6;
+	tempArray[7] = mp7;
+	tempArray[8] = NULL;
+
+	while (tempArray[n] != NULL) {
+		if (mp_init(pool, tempArray[n]) != MP_OKAY) {
+			res = MP_MEM;
+			break;
+		}
+		n++;
+	}
+
+	if (res == MP_MEM) {
+		n = 0;
+		while (tempArray[n] != NULL) {
+			mp_clear(tempArray[n]);
+			n++;
+		}
+	}
+	return res;		/* Assumed ok, if error flagged above. */
+}
+/******************************************************************************/
+/*
+	Reads a unsigned char array, assumes the msb is stored first [big endian]
+ */
+int32 mp_read_unsigned_bin (mp_int * a, unsigned char *b, int32 c)
+{
+	int32		res;
+
+/*
+	Make sure there are at least two digits.
+ */
+	if (a->alloc < 2) {
+		if ((res = mp_grow(a, 2)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+/*
+	Zero the int32.
+ */
+	mp_zero (a);
+
+/*
+	read the bytes in
+ */
+	while (c-- > 0) {
+		if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) {
+			return res;
+		}
+
+#ifndef MP_8BIT
+		a->dp[0] |= *b++;
+		a->used += 1;
+#else
+		a->dp[0] = (*b & MP_MASK);
+		a->dp[1] |= ((*b++ >> 7U) & 1);
+		a->used += 2;
+#endif /* MP_8BIT */
+	}
+	mp_clamp (a);
+	return MP_OKAY;
+}
+
+/******************************************************************************/
+/* 
+	Compare two ints (signed)
+ */
+int32 mp_cmp (mp_int * a, mp_int * b)
+{
+/*
+	compare based on sign
+ */
+	if (a->sign != b->sign) {
+		if (a->sign == MP_NEG) {
+			return MP_LT;
+		} else {
+			return MP_GT;
+		}
+	}
+
+/*
+	compare digits
+ */
+	if (a->sign == MP_NEG) {
+		/* if negative compare opposite direction */
+		return mp_cmp_mag(b, a);
+	} else {
+		return mp_cmp_mag(a, b);
+	}
+}
+
+/******************************************************************************/
+/*
+	Store in unsigned [big endian] format.
+*/
+int32 mp_to_unsigned_bin(psPool_t *pool, mp_int * a, unsigned char *b)
+{
+	int32			x, res;
+	mp_int		t;
+
+	if ((res = mp_init_copy(pool, &t, a)) != MP_OKAY) {
+		return res;
+	}
+
+	x = 0;
+	while (mp_iszero (&t) == 0) {
+#ifndef MP_8BIT
+		b[x++] = (unsigned char) (t.dp[0] & 255);
+#else
+		b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
+#endif /* MP_8BIT */
+		if ((res = mp_div_2d (pool, &t, 8, &t, NULL)) != MP_OKAY) {
+			mp_clear (&t);
+			return res;
+		}
+	}
+	bn_reverse (b, x);
+	mp_clear (&t);
+	return MP_OKAY;
+}
+
+void _mp_clear_multi(mp_int *mp0, mp_int *mp1, mp_int *mp2, mp_int *mp3,
+				  mp_int *mp4, mp_int *mp5, mp_int *mp6, mp_int *mp7)
+{
+	int32		n		= 0;		/* Number of ok inits */
+
+	mp_int	*tempArray[9];
+	
+	tempArray[0] = mp0;
+	tempArray[1] = mp1;
+	tempArray[2] = mp2;
+	tempArray[3] = mp3;
+	tempArray[4] = mp4;
+	tempArray[5] = mp5;
+	tempArray[6] = mp6;
+	tempArray[7] = mp7;
+	tempArray[8] = NULL;
+
+	for (n = 0; tempArray[n] != NULL; n++) {
+		mp_clear(tempArray[n]);
+	}
+}
+
+/******************************************************************************/
+/*
+	Init a new mp_int.
+*/
+int32 mp_init (psPool_t *pool, mp_int * a)
+{
+	int32		i;
+/*
+	allocate memory required and clear it
+ */
+	a->dp = OPT_CAST(mp_digit) psMalloc(pool, sizeof (mp_digit) * MP_PREC);
+	if (a->dp == NULL) {
+		return MP_MEM;
+	}
+
+/*
+	set the digits to zero
+ */
+	for (i = 0; i < MP_PREC; i++) {
+		a->dp[i] = 0;
+	}
+/*
+	set the used to zero, allocated digits to the default precision and sign
+	to positive
+ */
+	a->used  = 0;
+	a->alloc = MP_PREC;
+	a->sign  = MP_ZPOS;
+
+	return MP_OKAY;
+}
+
+/******************************************************************************/
+/*
+	clear one (frees).
+ */
+void mp_clear (mp_int * a)
+{
+	int32		i;
+/*
+	only do anything if a hasn't been freed previously
+ */
+	if (a->dp != NULL) {
+/*
+		first zero the digits
+ */
+		for (i = 0; i < a->used; i++) {
+			a->dp[i] = 0;
+		}
+
+		/* free ram */
+		psFree (a->dp);
+
+/*
+		reset members to make debugging easier
+ */
+		a->dp		= NULL;
+		a->alloc	= a->used = 0;
+		a->sign		= MP_ZPOS;
+	}
+}
+
+/******************************************************************************/
+/*
+	Get the size for an unsigned equivalent.
+ */
+int32 mp_unsigned_bin_size (mp_int * a)
+{
+	int32		size = mp_count_bits (a);
+
+	return	(size / 8 + ((size & 7) != 0 ? 1 : 0));
+}
+
+/******************************************************************************/
+/*
+	Trim unused digits 
+
+	This is used to ensure that leading zero digits are trimed and the 
+	leading "used" digit will be non-zero. Typically very fast.  Also fixes 
+	the sign if there are no more leading digits
+*/
+void mp_clamp (mp_int * a)
+{
+/*
+	decrease used while the most significant digit is zero.
+ */
+	while (a->used > 0 && a->dp[a->used - 1] == 0) {
+		--(a->used);
+	}
+
+/*
+	reset the sign flag if used == 0
+ */
+	if (a->used == 0) {
+		a->sign = MP_ZPOS;
+	}
+}
+
+/******************************************************************************/
+/*
+	Shift left by a certain bit count.
+ */
+int32 mp_mul_2d (mp_int * a, int32 b, mp_int * c)
+{
+	mp_digit	d;
+	int32			res;
+
+/*
+	Copy
+ */
+	if (a != c) {
+		if ((res = mp_copy (a, c)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+	if (c->alloc < (int32)(c->used + b/DIGIT_BIT + 1)) {
+		if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+/*
+	Shift by as many digits in the bit count
+ */
+	if (b >= (int32)DIGIT_BIT) {
+		if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+/*
+	shift any bit count < DIGIT_BIT
+ */
+	d = (mp_digit) (b % DIGIT_BIT);
+	if (d != 0) {
+		register mp_digit *tmpc, shift, mask, r, rr;
+		register int32 x;
+
+/*
+		bitmask for carries
+ */
+		mask = (((mp_digit)1) << d) - 1;
+
+/*
+		shift for msbs
+ */
+		shift = DIGIT_BIT - d;
+
+		/* alias */
+		tmpc = c->dp;
+
+		/* carry */
+		r = 0;
+		for (x = 0; x < c->used; x++) {
+/*
+			get the higher bits of the current word
+ */
+			rr = (*tmpc >> shift) & mask;
+
+/*
+			shift the current word and OR in the carry
+ */
+			*tmpc = ((*tmpc << d) | r) & MP_MASK;
+			++tmpc;
+
+/*
+			set the carry to the carry bits of the current word
+ */
+			r = rr;
+		}
+
+/*
+		set final carry
+ */
+		if (r != 0) {
+			c->dp[(c->used)++] = r;
+		}
+	}
+	mp_clamp (c);
+	return MP_OKAY;
+}
+
+/******************************************************************************/
+/* 
+	Set to zero.
+ */
+void mp_zero (mp_int * a)
+{
+	int			n;
+	mp_digit	*tmp;
+
+	a->sign = MP_ZPOS;
+	a->used = 0;
+
+	tmp = a->dp;
+	for (n = 0; n < a->alloc; n++) {
+		*tmp++ = 0;
+	}
+}
+
+#ifdef MP_LOW_MEM
+#define TAB_SIZE 32
+#else
+#define TAB_SIZE 256
+#endif /* MP_LOW_MEM */
+
+/******************************************************************************/
+/*
+	Compare maginitude of two ints (unsigned).
+ */
+int32 mp_cmp_mag (mp_int * a, mp_int * b)
+{
+	int32			n;
+	mp_digit	*tmpa, *tmpb;
+
+/*
+	compare based on # of non-zero digits
+ */
+	if (a->used > b->used) {
+		return MP_GT;
+	}
+
+	if (a->used < b->used) {
+		return MP_LT;
+	}
+
+	/* alias for a */
+	tmpa = a->dp + (a->used - 1);
+
+	/* alias for b */
+	tmpb = b->dp + (a->used - 1);
+
+/*
+	compare based on digits
+ */
+	for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
+		if (*tmpa > *tmpb) {
+			return MP_GT;
+		}
+
+		if (*tmpa < *tmpb) {
+			return MP_LT;
+		}
+	}
+	return MP_EQ;
+}
+
+/******************************************************************************/
+/*
+	computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
+
+	Uses a left-to-right k-ary sliding window to compute the modular 
+	exponentiation. The value of k changes based on the size of the exponent.
+
+	Uses Montgomery or Diminished Radix reduction [whichever appropriate]
+*/
+int32 mp_exptmod(psPool_t *pool, mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
+{
+
+/*
+	modulus P must be positive
+ */
+	if (P->sign == MP_NEG) {
+		return MP_VAL;
+	}
+
+/*
+	if exponent X is negative we have to recurse
+ */
+	if (X->sign == MP_NEG) {
+		mp_int tmpG, tmpX;
+		int32 err;
+
+/*
+		first compute 1/G mod P
+ */
+		if ((err = mp_init(pool, &tmpG)) != MP_OKAY) {
+			return err;
+		}
+		if ((err = mp_invmod(pool, G, P, &tmpG)) != MP_OKAY) {
+			mp_clear(&tmpG);
+			return err;
+		}
+
+/*
+		now get |X|
+ */
+		if ((err = mp_init(pool, &tmpX)) != MP_OKAY) {
+			mp_clear(&tmpG);
+			return err;
+		}
+		if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
+			mp_clear(&tmpG);
+			mp_clear(&tmpX);
+			return err;
+		}
+
+/*
+		and now compute (1/G)**|X| instead of G**X [X < 0]
+ */
+		err = mp_exptmod(pool, &tmpG, &tmpX, P, Y);
+		mp_clear(&tmpG);
+		mp_clear(&tmpX);
+		return err;
+	}
+
+/*
+	if the modulus is odd or dr != 0 use the fast method
+ */
+	if (mp_isodd (P) == 1) {
+		return mp_exptmod_fast (pool, G, X, P, Y, 0);
+	} else {
+/*
+		no exptmod for evens
+ */
+		return MP_VAL;
+	}
+}
+
+/******************************************************************************/
+/*
+	Call only from mp_exptmod to make sure this fast version qualifies
+*/
+static int32 mp_exptmod_fast(psPool_t *pool, mp_int * G, mp_int * X,
+				mp_int * P, mp_int * Y, int32 redmode)
+{
+	mp_int		M[TAB_SIZE], res;
+	mp_digit	buf, mp;
+	int32		err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
+
+
+/*
+	use a pointer to the reduction algorithm.  This allows us to use
+	one of many reduction algorithms without modding the guts of
+	the code with if statements everywhere.
+ */
+	int32	(*redux)(mp_int*,mp_int*,mp_digit);
+
+/*
+	find window size
+ */
+	x = mp_count_bits (X);
+	if (x <= 7) {
+		winsize = 2;
+	} else if (x <= 36) {
+		winsize = 3;
+	} else if (x <= 140) {
+		winsize = 4;
+	} else if (x <= 450) {
+		winsize = 5;
+	} else if (x <= 1303) {
+		winsize = 6;
+	} else if (x <= 3529) {
+		winsize = 7;
+	} else {
+		winsize = 8;
+	}
+
+#ifdef MP_LOW_MEM
+	if (winsize > 5) {
+		winsize = 5;
+	}
+#endif
+
+/*
+	init M array
+	init first cell
+ */
+	if ((err = mp_init(pool, &M[1])) != MP_OKAY) {
+		return err;
+	}
+
+/*
+	now init the second half of the array
+ */
+	for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
+		if ((err = mp_init(pool, &M[x])) != MP_OKAY) {
+			for (y = 1<<(winsize-1); y < x; y++) {
+				mp_clear(&M[y]);
+			}
+			mp_clear(&M[1]);
+			return err;
+		}
+	}
+
+
+/*
+	now setup montgomery
+ */
+	if ((err = mp_montgomery_setup(P, &mp)) != MP_OKAY) {
+		goto LBL_M;
+	}
+
+/*
+	automatically pick the comba one if available
+ */
+	if (((P->used * 2 + 1) < MP_WARRAY) &&
+			P->used < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
+		redux = fast_mp_montgomery_reduce;
+	} else {
+/*
+		use slower baseline Montgomery method
+ */
+		redux = mp_montgomery_reduce;
+	}
+
+/*
+	setup result
+ */
+	if ((err = mp_init(pool, &res)) != MP_OKAY) {
+		goto LBL_M;
+	}
+
+/*
+	create M table. The first half of the table is not computed
+	though accept for M[0] and M[1]
+*/
+
+/*
+	now we need R mod m
+ */
+	if ((err = mp_montgomery_calc_normalization(&res, P)) != MP_OKAY) {
+		goto LBL_RES;
+	}
+
+/*
+	now set M[1] to G * R mod m
+ */
+	if ((err = mp_mulmod(pool, G, &res, P, &M[1])) != MP_OKAY) {
+		goto LBL_RES;
+	}
+
+/*
+	compute the value at M[1<<(winsize-1)] by squaring
+	M[1] (winsize-1) times
+*/
+	if ((err = mp_copy(&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
+		goto LBL_RES;
+	}
+
+	for (x = 0; x < (winsize - 1); x++) {
+		if ((err = mp_sqr(pool, &M[1 << (winsize - 1)],
+				&M[1 << (winsize - 1)])) != MP_OKAY) {
+			goto LBL_RES;
+		}
+		if ((err = redux(&M[1 << (winsize - 1)], P, mp)) != MP_OKAY) {
+			goto LBL_RES;
+		}
+	}
+
+/*
+	create upper table
+ */
+	for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
+		if ((err = mp_mul(pool, &M[x - 1], &M[1], &M[x])) != MP_OKAY) {
+			goto LBL_RES;
+		}
+		if ((err = redux(&M[x], P, mp)) != MP_OKAY) {
+			goto LBL_RES;
+		}
+	}
+
+/*
+	set initial mode and bit cnt
+ */
+	mode   = 0;
+	bitcnt = 1;
+	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;
+	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);
+		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
+ */
+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++) {
+/*
+		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;
+}
+
+/******************************************************************************/
+/*
+	FUTURE - this is only needed by the SSH code, SLOW or not, because RSA
+	exponents are always odd.
+*/
+int32 mp_invmodSSH(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;
+}
+
+/******************************************************************************/
+
+/*
+ *	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) != MP_EQ) {
+		res = MP_VAL;
+		goto LBL_ERR;
+	}
+
+/*
+	b is now the inverse
+ */
+	neg = a->sign;
+	while (D.sign == MP_NEG) {
+		if ((res = mp_add(&D, b, &D)) != MP_OKAY) {
+		goto LBL_ERR;
+		}
+	}
+	mp_exch(&D, c);
+	c->sign = neg;
+	res = MP_OKAY;
+
+LBL_ERR:_mp_clear_multi(&x, &y, &u, &v, &B, &D, NULL, NULL);
+	return res;
+}
+
+/******************************************************************************/
+/*
+	d = a + b (mod c)
+ */
+int32 mp_addmod (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_add (a, b, &t)) != MP_OKAY) {
+		mp_clear (&t);
+		return res;
+	}
+	res = mp_mod (pool, &t, c, d);
+	mp_clear (&t);
+	return res;
+}
+
+/******************************************************************************/
+/*
+	shrink a bignum
+ */
+int32 mp_shrink (mp_int * a)
+{
+	mp_digit *tmp;
+
+	if (a->alloc != a->used && a->used > 0) {
+		if ((tmp = psRealloc(a->dp, sizeof (mp_digit) * a->used)) == NULL) {
+			return MP_MEM;
+		}
+		a->dp    = tmp;
+		a->alloc = a->used;
+	}
+	return MP_OKAY;
+}
+
+/* single digit subtraction */
+int32 mp_sub_d (mp_int * a, mp_digit b, mp_int * c)
+{
+	mp_digit *tmpa, *tmpc, mu;
+	int32       res, ix, oldused;
+
+	/* grow c as required */
+	if (c->alloc < a->used + 1) {
+		if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+	/* if a is negative just do an unsigned
+	* addition [with fudged signs]
+	*/
+	if (a->sign == MP_NEG) {
+		a->sign = MP_ZPOS;
+		res     = mp_add_d(a, b, c);
+		a->sign = c->sign = MP_NEG;
+		return res;
+	}
+
+	/* setup regs */
+	oldused = c->used;
+	tmpa    = a->dp;
+	tmpc    = c->dp;
+
+	/* if a <= b simply fix the single digit */
+	if ((a->used == 1 && a->dp[0] <= b) || a->used == 0) {
+		if (a->used == 1) {
+			*tmpc++ = b - *tmpa;
+		} else {
+			*tmpc++ = b;
+		}
+		ix      = 1;
+
+		/* negative/1digit */
+		c->sign = MP_NEG;
+		c->used = 1;
+	} else {
+		/* positive/size */
+		c->sign = MP_ZPOS;
+		c->used = a->used;
+
+		/* subtract first digit */
+		*tmpc    = *tmpa++ - b;
+		mu       = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
+		*tmpc++ &= MP_MASK;
+
+		/* handle rest of the digits */
+		for (ix = 1; ix < a->used; ix++) {
+			*tmpc    = *tmpa++ - mu;
+			mu       = *tmpc >> (sizeof(mp_digit) * CHAR_BIT - 1);
+			*tmpc++ &= MP_MASK;
+		}
+	}
+
+	/* zero excess digits */
+	while (ix++ < oldused) {
+		*tmpc++ = 0;
+	}
+	mp_clamp(c);
+	return MP_OKAY;
+}
+
+/* single digit addition */
+int32 mp_add_d (mp_int * a, mp_digit b, mp_int * c)
+{
+	int32     res, ix, oldused;
+	mp_digit *tmpa, *tmpc, mu;
+
+	/* grow c as required */
+	if (c->alloc < a->used + 1) {
+		if ((res = mp_grow(c, a->used + 1)) != MP_OKAY) {
+			return res;
+		}
+	}
+
+	/* if a is negative and |a| >= b, call c = |a| - b */
+	if (a->sign == MP_NEG && (a->used > 1 || a->dp[0] >= b)) {
+		/* temporarily fix sign of a */
+		a->sign = MP_ZPOS;
+
+		/* c = |a| - b */
+		res = mp_sub_d(a, b, c);
+
+		/* fix sign  */
+		a->sign = c->sign = MP_NEG;
+		return res;
+	}
+
+	/* old number of used digits in c */
+	oldused = c->used;
+
+	/* sign always positive */
+	c->sign = MP_ZPOS;
+
+	/* source alias */
+	tmpa    = a->dp;
+
+	/* destination alias */
+	tmpc    = c->dp;
+
+	/* if a is positive */
+	if (a->sign == MP_ZPOS) {
+		/* add digit, after this we're propagating the carry */
+		*tmpc   = *tmpa++ + b;
+		mu      = *tmpc >> DIGIT_BIT;
+		*tmpc++ &= MP_MASK;
+
+		/* now handle rest of the digits */
+		for (ix = 1; ix < a->used; ix++) {
+			*tmpc   = *tmpa++ + mu;
+			mu      = *tmpc >> DIGIT_BIT;
+			*tmpc++ &= MP_MASK;
+		}
+		/* set final carry */
+		ix++;
+		*tmpc++  = mu;
+
+		/* setup size */
+		c->used = a->used + 1;
+	} else {
+		/* a was negative and |a| < b */
+		c->used  = 1;
+
+		/* the result is a single digit */
+		if (a->used == 1) {
+			*tmpc++  =  b - a->dp[0];
+		} else {
+			*tmpc++  =  b;
+		}
+
+		/* setup count so the clearing of oldused
+		* can fall through correctly
+		*/
+		ix       = 1;
+	}
+
+	/* now zero to oldused */
+	while (ix++ < oldused) {
+		*tmpc++ = 0;
+	}
+	mp_clamp(c);
+	return MP_OKAY;
+}
+
+
+/******************************************************************************/
+
+#endif /* USE_MPI2 */