vmbigintegersclass.cpp

TOOL (Tiny Object Oriented Language) is an easily-embedded, object-oriented, C++-like-language inter

	p1 = qhat * v1;
	t = p0 + VM_TOHIGH(VM_LOHALF(p1));
	uu[0] -= t;
	if (uu[0] > VM_MAX_DIGIT - t)
		uu[1]--;	/* Borrow */
	uu[1] -= VM_HIHALF(p1);
}



VM_DIGIT_T VMBigIntegers::mpSubtract(VM_DIGIT_T w[], VM_DIGIT_T u[], VM_DIGIT_T v[], 
			   unsigned int ndigits)
{
	/*	Calculates w = u - v where u >= v
		w, u, v are multiprecision integers of ndigits each
		Returns 0 if OK, or 1 if v > u.

		Ref: Knuth Vol 2 Ch 4.3.1 p 267 Algorithm S.
	*/

	VM_DIGIT_T k;
	unsigned int j;

	/* Step S1. Initialise */
	k = 0;

	for (j = 0; j < ndigits; j++)
	{
		/*	Step S2. Subtract digits w_j = (u_j - v_k - k)
			Set k = 1 if borrow occurs.
		*/
		w[j] = u[j] - k;
		if (w[j] > VM_MAX_DIGIT - k)
			k = 1;
		else
			k = 0;
		
		w[j] -= v[j];
		if (w[j] > VM_MAX_DIGIT - v[j])
			k++;

	}	/* Step S3. Loop on j */

	return k;	/* Should be zero if u >= v */
}


unsigned int VMBigIntegers::mpSizeof(VM_DIGIT_T a[], unsigned int ndigits)
{	/* Returns size of significant digits in a */
	
	while(ndigits--)
	{
		if (a[ndigits] != 0)
			return (++ndigits);
	}
	return 0;
}


VM_DIGIT_T VMBigIntegers::mpShortSub(VM_DIGIT_T w[], VM_DIGIT_T u[], VM_DIGIT_T v, 
			   unsigned int ndigits)
{
	/*	Calculates w = u - v
		where w, u are multiprecision integers of ndigits each
		and v is a single precision digit.
		Returns borrow: 0 if u >= v, or 1 if v > u.

		Ref: Derived from Knuth Algorithm S.
	*/

	VM_DIGIT_T k;
	unsigned int j;

	k = 0;

	/* Subtract v from first digit of u */
	w[0] = u[0] - v;
	if (w[0] > VM_MAX_DIGIT - v)
		k = 1;
	else
		k = 0;

	/* Subtract borrow from subsequent digits */
	for (j = 1; j < ndigits; j++)
	{
		w[j] = u[j] - k;
		if (w[j] > VM_MAX_DIGIT - k)
			k = 1;
		else
			k = 0;
	}

	return k;
}


VM_DIGIT_T VMBigIntegers::mpShortMult(VM_DIGIT_T w[], VM_DIGIT_T u[], VM_DIGIT_T v, 
					unsigned int ndigits)
{
	/*	Computes product w = u * v
		Returns overflow k
		where w, u are multiprecision integers of ndigits each
		and v, k are single precision digits

		Ref: Knuth Algorithm M.
	*/

	VM_DIGIT_T k, t[2];
	unsigned int j;

	/*
	for (j = 0; j < ndigits; j++)
		w[j] = 0;
	*/

	if (v == 0) return 0;

	k = 0;
	for (j = 0; j < ndigits; j++)
	{
		/* t = x_i * v */
		spMultiply(t, u[j], v);
		/* w_i = VM_LOHALF(t) + carry */
		w[j] = t[0] + k;
		/* Overflow? */
		if (w[j] < k)
			t[1]++;
		/* Carry forward VM_HIHALF(t) */
		k = t[1];
	}

	return k;
}


VM_DIGIT_T VMBigIntegers::mpShortMod(VM_DIGIT_T a[], VM_DIGIT_T d, unsigned int ndigits)
{
	/*	Calculates r = a mod d
		where a is a multiprecision integer of ndigits
		and r, d are single precision digits
		using bit-by-bit method from left to right.

		Use remainder from divide function.
	*/

	VM_DIGIT_T q[VM_MAX_DIG_LEN * 2];
	VM_DIGIT_T r = 0;

	r= mpShortDiv(q, a, d, ndigits);

	mpSetZero(q, ndigits);

	return r;
}


VM_DIGIT_T VMBigIntegers::mpShortDiv(VM_DIGIT_T q[], VM_DIGIT_T u[], VM_DIGIT_T v, 
				   unsigned int ndigits)
{
	/*	Calculates quotient q = u div v
		Returns remainder r = u mod v
		where q, u are multiprecision integers of ndigits each
		and d, v are single precision digits.

		Makes no assumptions about normalisation.
		
		Ref: Knuth Vol 2 Ch 4.3.1 Exercise 16 p625
	*/
	unsigned int j;
	VM_DIGIT_T t[2], r;
	unsigned int shift;
	VM_DIGIT_T bitmask, overflow, *uu;

	if (ndigits == 0) return 0;
	if (v == 0)	return 0;	/* Divide by zero error */

	/*	Normalise first */
	/*	Requires high bit of V
		to be set, so find most signif. bit then shift left,
		i.e. d = 2^shift, u' = u * d, v' = v * d.
	*/
	bitmask = VM_HIBITMASK;
	for (shift = 0; shift < VM_BITS_PER_DIGIT; shift++)
	{
		if (v & bitmask)
			break;
		bitmask >>= 1;
	}

	v <<= shift;
	overflow = mpShiftLeft(q, u, shift, ndigits);
	uu = q;
	
	/* Step S1 - modified for extra digit. */
	r = overflow;	/* New digit Un */
	j = ndigits;
	while (j--)
	{
		/* Step S2. */
		t[1] = r;
		t[0] = uu[j];
		overflow = spDivide(&q[j], &r, t, v);
	}

	/* Unnormalise */
	r >>= shift;
	
	return r;
}


int VMBigIntegers::mpShortCmp(VM_DIGIT_T a[], VM_DIGIT_T b, unsigned int ndigits)
{
	/* Returns sign of (a - b) where b is a single digit */

	unsigned int i;

	if (ndigits == 0) return 0;

	for (i = 1; i < ndigits; i++)
	{
		if (a[i] != 0)
			return 1;	/* GT */
	}

	if (a[0] < b)
		return -1;		/* LT */
	else if (a[0] > b)
		return 1;		/* GT */

	return 0;	/* EQ */
}


VM_DIGIT_T VMBigIntegers::mpShortAdd(VM_DIGIT_T w[], VM_DIGIT_T u[], VM_DIGIT_T v, 
			   unsigned int ndigits)
{
	/*	Calculates w = u + v
		where w, u are multiprecision integers of ndigits each
		and v is a single precision digit.
		Returns carry if overflow.

		Ref: Derived from Knuth Algorithm A.
	*/

	VM_DIGIT_T k;
	unsigned int j;

	k = 0;

	/* Add v to first digit of u */
	w[0] = u[0] + v;
	if (w[0] < v)
		k = 1;
	else
		k = 0;

	/* Add carry to subsequent digits */
	for (j = 1; j < ndigits; j++)
	{
		w[j] = u[j] + k;
		if (w[j] < k)
			k = 1;
		else
			k = 0;
	}

	return k;
}


VM_DIGIT_T VMBigIntegers::mpShiftRight(VM_DIGIT_T a[], VM_DIGIT_T b[],
	unsigned int x, unsigned int ndigits)
{	/* Computes a = b >> x */
	unsigned int i, y;
	VM_DIGIT_T mask, carry, nextcarry;

	/* Check input  - NB unspecified result */
	if (x >= VM_BITS_PER_DIGIT)
		return 0;

	/* Construct mask */
	mask = 0x1;
	for (i = 1; i < x; i++)
	{
		mask = (mask << 1) | mask;
	}
	if (x == 0) mask = 0x0;
	
	y = VM_BITS_PER_DIGIT - x;
	carry = 0;
	i = ndigits;
	while (i--)
	{
		nextcarry = (b[i] & mask) << y;
		a[i] = b[i] >> x | carry;
		carry = nextcarry;
	}

	return carry;
}


VM_DIGIT_T VMBigIntegers::mpShiftLeft(VM_DIGIT_T a[], VM_DIGIT_T b[],
	unsigned int x, unsigned int ndigits)
{	/* Computes a = b << x */
	unsigned int i, y;
	VM_DIGIT_T mask, carry, nextcarry;

	/* Check input - NB unspecified result */
	if (x >= VM_BITS_PER_DIGIT)
		return 0;

	/* Construct mask */
	mask = VM_HIBITMASK;
	for (i = 1; i < x; i++)
	{
		mask = (mask >> 1) | mask;
	}
	if (x == 0) mask = 0x0;
	
	y = VM_BITS_PER_DIGIT - x;
	carry = 0;
	for (i = 0; i < ndigits; i++)
	{
		nextcarry = (b[i] & mask) >> y;
		a[i] = b[i] << x | carry;
		carry = nextcarry;
	}

	return carry;
}


void VMBigIntegers::mpSetZero(VM_DIGIT_T a[], unsigned int ndigits)
{	/* Sets a = 0 */
	unsigned int i;
	
	for (i = 0; i < ndigits; i++)
	{
		a[i] = 0;
	}
}


void VMBigIntegers::mpSetEqual(VM_DIGIT_T a[], VM_DIGIT_T b[], unsigned int ndigits)
{	/* Sets a = b */
	unsigned int i;
	
	for (i = 0; i < ndigits; i++)
	{
		a[i] = b[i];
	}
}


void VMBigIntegers::mpSetDigit(VM_DIGIT_T a[], VM_DIGIT_T d, unsigned int ndigits)
{	/* Sets a = d where d is a single digit */
	unsigned int i;
	
	for (i = 1; i < ndigits; i++)
	{
		a[i] = 0;
	}
	a[0] = d;
}


int VMBigIntegers::mpMultiply(VM_DIGIT_T w[], VM_DIGIT_T u[], VM_DIGIT_T v[], 
					unsigned int ndigits)
{
	/*	Computes product w = u * v
		where u, v are multiprecision integers of ndigits each
		and w is a multiprecision integer of 2*ndigits

		Ref: Knuth Vol 2 Ch 4.3.1 p 268 Algorithm M.
	*/

	VM_DIGIT_T k, t[2];
	unsigned int i, j, m, n;

	m = n = ndigits;

	/* Step M1. Initialise */
	for (i = 0; i < 2 * m; i++)
		w[i] = 0;

	for (j = 0; j < n; j++)
	{
		/* Step M2. Zero multiplier? */
		if (v[j] == 0)
		{
			w[j + m] = 0;
		}
		else
		{
			/* Step M3. Initialise i */
			k = 0;
			for (i = 0; i < m; i++)
			{
				/* Step M4. Multiply and add */
				/* t = u_i * v_j + w_(i+j) + k */
				spMultiply(t, u[i], v[j]);

				t[0] += k;
				if (t[0] < k)
					t[1]++;
				t[0] += w[i+j];
				if (t[0] < w[i+j])
					t[1]++;

				w[i+j] = t[0];
				k = t[1];
			}	
			/* Step M5. Loop on i, set w_(j+m) = k */
			w[j+m] = k;
		}
	}	/* Step M6. Loop on j */

	return 0;
}


int VMBigIntegers::mpModulo(VM_DIGIT_T r[], VM_DIGIT_T u[], unsigned int udigits, 
			 VM_DIGIT_T v[], unsigned int vdigits)
{
	/*	Calculates r = u mod v
		where r, v are multiprecision integers of length vdigits
		and u is a multiprecision integer of length udigits.
		r may overlap v.

		Note that r here is only vdigits long, 
		whereas in mpDivide it is udigits long.

		Use remainder from mpDivide function.
	*/

	/* Double-length temp variable for divide fn */
	VM_DIGIT_T qq[VM_MAX_DIG_LEN * 2];
	/* Use a double-length temp for r to allow overlap of r and v */
	VM_DIGIT_T rr[VM_MAX_DIG_LEN * 2];

	/* rr[2n] = u[2n] mod v[n] */
	mpDivide(qq, rr, u, udigits, v, vdigits);

	mpSetEqual(r, rr, vdigits);
	mpSetZero(rr, udigits);
	mpSetZero(qq, udigits);

	return 0;
}


int VMBigIntegers::mpModMult(VM_DIGIT_T a[], VM_DIGIT_T x[], VM_DIGIT_T y[], 
			  VM_DIGIT_T m[], unsigned int ndigits)
{	/*	Computes a = (x * y) mod m */
	
	/* Double-length temp variable */
	VM_DIGIT_T p[VM_MAX_DIG_LEN * 2];

	/* Calc p[2n] = x * y */
	mpMultiply(p, x, y, ndigits);

	/* Then modulo */
	mpModulo(a, p, ndigits * 2, m, ndigits);

	mpSetZero(p, ndigits * 2);
	return 0;
}


int VMBigIntegers::mpModInv(VM_DIGIT_T inv[], VM_DIGIT_T u[], VM_DIGIT_T v[], unsigned int ndigits)
{	/*	Computes inv = u^(-1) mod v */
	/*	Ref: Knuth Algorithm X Vol 2 p 342
		ignoring u2, v2, t2
		and avoiding negative numbers.
	*/
	/* Allocate temp variables */
	VM_DIGIT_T u1[VM_MAX_DIG_LEN], u3[VM_MAX_DIG_LEN], v1[VM_MAX_DIG_LEN], v3[VM_MAX_DIG_LEN];
	VM_DIGIT_T t1[VM_MAX_DIG_LEN], t3[VM_MAX_DIG_LEN], q[VM_MAX_DIG_LEN];
	VM_DIGIT_T w[2*VM_MAX_DIG_LEN];
	/* TODO: CHECK LENGTHS HERE */
	int bIterations;

	/* Step X1. Initialise */
	mpSetDigit(u1, 1, ndigits);		/* u1 = 1 */
	mpSetEqual(u3, u, ndigits);		/* u3 = u */
	mpSetZero(v1, ndigits);			/* v1 = 0 */
	mpSetEqual(v3, v, ndigits);		/* v3 = v */

	bIterations = 1;	/* Remember odd/even iterations */
	while (!mpIsZero(v3, ndigits))		/* Step X2. Loop while v3 != 0 */
	{					/* Step X3. Divide and "Subtract" */
		mpDivide(q, t3, u3, ndigits, v3, ndigits);
						/* q = u3 / v3, t3 = u3 % v3 */
		mpMultiply(w, q, v1, ndigits);	/* w = q * v1 */
		mpAdd(t1, u1, w, ndigits);		/* t1 = u1 + w */

		/* Swap u1 = v1; v1 = t1; u3 = v3; v3 = t3 */
		mpSetEqual(u1, v1, ndigits);
		mpSetEqual(v1, t1, ndigits);
		mpSetEqual(u3, v3, ndigits);
		mpSetEqual(v3, t3, ndigits);

		bIterations = -bIterations;
	}

	if (bIterations < 0)
		mpSubtract(inv, v, u1, ndigits);	/* inv = v - u1 */
	else
		mpSetEqual(inv, u1, ndigits);	/* inv = u1 */

	/* Clear up */