vmbigintegersclass.cpp

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

	mpSetZero(u1, ndigits);
	mpSetZero(v1, ndigits);
	mpSetZero(t1, ndigits);
	mpSetZero(u3, ndigits);
	mpSetZero(v3, ndigits);
	mpSetZero(t3, ndigits);
	mpSetZero(q, ndigits);
	mpSetZero(w, 2*ndigits);

	return 0;
}


int VMBigIntegers::mpModExp(VM_DIGIT_T y[], VM_DIGIT_T x[], 
			VM_DIGIT_T e[], VM_DIGIT_T m[], unsigned int ndigits)
{	/*	Computes y = x^e mod m */
	/*	Binary left-to-right method
	*/
	VM_DIGIT_T mask;
	unsigned int n;
	
	if (ndigits == 0) return -1;

	/* Find second-most significant bit in e */
	n = mpSizeof(e, ndigits);
	for (mask = VM_HIBITMASK; mask > 0; mask >>= 1)
	{
		if (e[n-1] & mask)
			break;
	}
	VM_mpNEXTBITMASK(mask, n);

	/* Set y = x */
	mpSetEqual(y, x, ndigits);

	/* For bit j = k-2 downto 0 step -1 */
	while (n)
	{
		mpModMult(y, y, y, m, ndigits);		/* Square */
		if (e[n-1] & mask)
			mpModMult(y, y, x, m, ndigits);	/* Multiply */
		
		/* Move to next bit */
		VM_mpNEXTBITMASK(mask, n);
	}

	return 0;
}


int VMBigIntegers::mpIsZero(VM_DIGIT_T a[], unsigned int ndigits)
{
	/*	Returns true if a == 0, else false
	*/

	unsigned int i;
	if (ndigits == 0) return -1;

	for (i = 0; i < ndigits; i++)	/* Start at lsb */
	{
		if (a[i] != 0)
			return 0;	/* False */
	}

	return (!0);	/* True */
}


int VMBigIntegers::mpIsPrime(VM_DIGIT_T w[], unsigned int ndigits, unsigned int t)
{	/*	Returns true if w is a probable prime 
		Carries out t iterations
		(Use t = 50 for DSS Standard) */
	/*	Uses Rabin-Miller Probabilistic Primality Test,
		Ref: FIPS-186-2 Appendix 2.
		Also Schneier 2nd ed p 260 & Knuth Vol 2, p 379.
	*/

	/* Temp big digits */
	VM_DIGIT_T m[VM_MAX_DIG_LEN], a[VM_MAX_DIG_LEN], b[VM_MAX_DIG_LEN];
	VM_DIGIT_T z[VM_MAX_DIG_LEN], w1[VM_MAX_DIG_LEN];
	VM_DIGIT_T j[VM_MAX_DIG_LEN];
	unsigned int i;
	int failed, isprime;

	/* Check the obvious */
	if (VM_mpISEVEN(w, ndigits))
		return 0;

	/*	First check for small primes */
	for (i = 0; i < VM_N_SMALL_PRIMES_2; i++)
	{
		if (mpShortMod(w ,VM_SMALL_PRIMES_2[i], ndigits) == 0)
			return 0;	/* Failed */
	}

	/*	Now do Rabin-Miller  */
	/*	Step 2. Find a and m where w = 1 + (2^a)m
		m is odd and 2^a is largest power of 2 dividing w - 1 */
	mpShortSub(w1, w, 1, ndigits);	/* Store w1 = w - 1 */
	mpSetEqual(m, w1, ndigits);		/* Set m = w - 1 */
	/* for (a = 0; iseven(m); a++) */
	for (mpSetZero(a, ndigits); VM_mpISEVEN(m, ndigits); 
		mpShortAdd(a, a, 1, ndigits))
	{	/* Divide by 2 until m is odd */
		mpShiftRight(m, m, 1, ndigits);
	}

	/* assert((1 << a) * m + 1 == w); */

	isprime = 1;
	for (i = 0; i < t; i++)
	{
		failed = 1;	/* Assume fail unless passed in loop */
		/* Step 3. Generate random integer 1 < b < w */
		mpSetZero(b, ndigits);
		do
		{
			b[0] = spPseudoRand(2, VM_MAX_DIGIT);
		} while (mpCompare(b, w, ndigits) >= 0);

		/* assert(1 < b && b < w); */

		/* Step 4. Set j = 0 and z = b^m mod w */
		mpSetZero(j, ndigits);
		mpModExp(z, b, m, w, ndigits);
		do
		{
			/* Step 5. If j = 0 and z = 1, or if z = w - 1 */
			/* i.e. if ((j == 0 && z == 1) || (z == w - 1)) */
			if ((mpIsZero(j, ndigits) 
				&& mpShortCmp(z, 1, ndigits) == 0)
				|| (mpCompare(z, w1, ndigits) == 0))
			{	/* Passes on this loop  - go to Step 9 */
				failed = 0;
				break;
			}

			/* Step 6. If j > 0 and z = 1 */
			if (!mpIsZero(j, ndigits) 
				&& (mpShortCmp(z, 1, ndigits) == 0))
			{	/* Fails - go to Step 8 */
				failed = 1;
				break;
			}

			/* Step 7. j = j + 1. If j < a set z = z^2 mod w */
			mpShortAdd(j, j, 1, ndigits);
			if (mpCompare(j, a, ndigits) < 0)
				mpModMult(z, z, z, w, ndigits);
			/* Loop: if j < a go to Step 5 */
		} while (mpCompare(j, a, ndigits) < 0);

		if (failed)
		{	/* Step 8. Not a prime - stop */
			isprime = 0;
			break;
		}
	}	/* Step 9. Go to Step 3 until i >= n */
	/* Else, if i = n, w is probably prime => success */

	/* Clean up */
	mpSetZero(m, ndigits);
	mpSetZero(a, ndigits);
	mpSetZero(b, ndigits);
	mpSetZero(z, ndigits);
	mpSetZero(w1, ndigits);

	return isprime;
}


VM_DIGIT_T VMBigIntegers::mpHalfMod(VM_DIGIT_T a[], VM_HALF_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.

		Ref: Derived from principles in PGP by Phil Zimmermann
		Note: This method will only work until r <<= 1 overflows.
		i.e. for d < VM_HIBITMASK, but we keep HALF_DIGIT
		limit for safety 
		(and also because we don't have a 31/32nds_digit).
	*/

	VM_DIGIT_T mask = VM_HIBITMASK;
	VM_DIGIT_T r = 0;

	if (ndigits == 0) return 0;

	while (ndigits)
	{	/* Work from left to right */

		r <<= 1;	/* Multiply remainder by 2 */

		/* Look at current bit */
		if (a[ndigits-1] & mask)
			r++;
		if (r >= d)
			r -= d;

		/* Move to next bit */
		if (mask == 1)
		{
			mask = VM_HIBITMASK;
			ndigits--;
		}
		else
			mask >>= 1;
	}
	return r;
}


VM_DIGIT_T VMBigIntegers::mpHalfDiv(VM_DIGIT_T q[], VM_DIGIT_T u[], VM_HALF_DIGIT_T v, 
				   unsigned int ndigits)
{
	return mpHalfDivK(q, u, v, ndigits);
}


VM_DIGIT_T VMBigIntegers::mpHalfDivK(VM_DIGIT_T q[], VM_DIGIT_T u[], VM_HALF_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
		
		d must be <= VM_MAX_HALF_DIGIT

		Ref: Knuth Vol 2 Ch 4.3.1 Exercise 16 p625
	*/
	unsigned int j;
	VM_DIGIT_T t, r, qHigh, qLow;

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

	/* Step S1. */
	r = 0;
	j = ndigits;
	while (j--)
	{
		/* Step S2. */
		t = VM_TOHIGH(r) | VM_HIHALF(u[j]);
		qHigh = t / v;
		r = t - qHigh * v;

		t = VM_TOHIGH(r) | VM_LOHALF(u[j]);
		qLow = t / v;
		r = t - qLow * v;

		q[j] = VM_TOHIGH(qHigh) | qLow;
	}

	return r;
}


VM_DIGIT_T VMBigIntegers::mpHalfDivZ(VM_DIGIT_T q[], VM_DIGIT_T a[], VM_HALF_DIGIT_T d, 
				   unsigned int ndigits)
{
	/*	Calculates quotient q = a div d
		Returns remainder r = a mod d
		where q, a are multiprecision integers of ndigits each
		and d, r are single precision digits
		using bit-by-bit method from left to right.

		d must be <= VM_MAX_HALF_DIGIT

		Ref: Principles in PGP by Phil Zimmermann
	*/

	VM_DIGIT_T mask = VM_HIBITMASK;
	VM_DIGIT_T r = 0;
	unsigned int i;

	if (ndigits == 0) return 0;

	/* Initialise quotient */
	for (i = 0; i < ndigits; i++)
		q[i] = 0;
	
	while (ndigits)
	{	/* Work from left to right */

		r <<= 1;	/* Multiply remainder by 2 */

		/* Look at current bit */
		if (a[ndigits-1] & mask)
			r++;
		if (r >= d)
		{
			r -= d;
			q[ndigits-1] |= mask;
		}

		/* Move to next bit */
		if (mask == 1)
		{
			mask = VM_HIBITMASK;
			ndigits--;
		}
		else
			mask >>= 1;
	}
	return r;
}


int VMBigIntegers::mpGcd(VM_DIGIT_T g[], VM_DIGIT_T x[], VM_DIGIT_T y[], unsigned int ndigits)
{	
	/* Computes g = gcd(x, y) */
	/* Ref: Schneier  */

	/*	NB This function requires temp storage 
	*/
	VM_DIGIT_T yy[VM_MAX_DIG_LEN], xx[VM_MAX_DIG_LEN];
	
	mpSetZero(yy, VM_MAX_DIG_LEN);
	mpSetZero(xx, VM_MAX_DIG_LEN);
	mpSetEqual(xx, x, ndigits);
	mpSetEqual(yy, y, ndigits);

	mpSetEqual(g, yy, ndigits);		/* g = y */
	
	while (!mpIsZero(xx, ndigits))	/* while (x > 0) */
	{
		mpSetEqual(g, xx, ndigits);	/* g = x */
		mpModulo(xx, yy, ndigits, xx, ndigits);	/* x = y mod x */
		mpSetEqual(yy, g, ndigits);	/* y = g; */
	}

	mpSetZero(xx, ndigits);
	mpSetZero(yy, ndigits);

	return 0;	/* gcd is in g */
}


int VMBigIntegers::mpEqual(VM_DIGIT_T a[], VM_DIGIT_T b[], unsigned int ndigits)
{
	/*	Returns true if a == b, else false
	*/

	if (ndigits == 0) return -1;

	while (ndigits--)
	{
		if (a[ndigits] != b[ndigits])
			return 0;	/* False */
	}

	return (!0);	/* True */
}


int VMBigIntegers::mpDivide(VM_DIGIT_T q[], VM_DIGIT_T r[], VM_DIGIT_T u[],
	unsigned int udigits, VM_DIGIT_T v[], unsigned int vdigits)
{	/*	Computes quotient q = u / v and remainder r = u mod v
		where q, r, u are multiple precision digits
		all of udigits and the divisor v is vdigits.

		Ref: Knuth Vol 2 Ch 4.3.1 p 272 Algorithm D.

		Do without extra storage space, i.e. use r[] for
		normalised u[], unnormalise v[] at end, and cope with
		extra digit Uj+n added to u after normalisation.

		WARNING: this trashes q and r first, so cannot do
		u = u / v or v = u mod v.
	*/
	unsigned int shift;
	int n, m, j;
	VM_DIGIT_T bitmask, overflow;
	VM_DIGIT_T qhat, rhat, t[2];
	VM_DIGIT_T *uu, *ww;
	int qhatOK, cmp;

	/* Clear q and r */
	mpSetZero(q, udigits);
	mpSetZero(r, udigits);

	/* Work out exact sizes of u and v */
	n = (int)mpSizeof(v, vdigits);
	m = (int)mpSizeof(u, udigits);
	m -= n;

	/* Catch special cases */
	if (n == 0)
		return -1;	/* Error: divide by zero */

	if (n == 1)
	{	/* Use short division instead */
		r[0] = mpShortDiv(q, u, v[0], udigits);
		return 0;
	}

	if (m < 0)
	{	/* v > u, so just set q = 0 and r = u */
		mpSetEqual(r, u, udigits);
		return 0;
	}

	if (m == 0)
	{	/* u and v are the same length */
		cmp = mpCompare(u, v, (unsigned int)n);
		if (cmp < 0)
		{	/* v > u, as above */
			mpSetEqual(r, u, udigits);
			return 0;
		}
		else if (cmp == 0)
		{	/* v == u, so set q = 1 and r = 0 */
			mpSetDigit(q, 1, udigits);
			return 0;
		}
	}

	/*	In Knuth notation, we have:
		Given
		u = (Um+n-1 ... U1U0)
		v = (Vn-1 ... V1V0)
		Compute
		q = u/v = (QmQm-1 ... Q0)
		r = u mod v = (Rn-1 ... R1R0)
	*/

	/*	Step D1. Normalise */
	/*	Requires high bit of Vn-1
		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[n-1] & bitmask)
			break;
		bitmask >>= 1;
	}

	/* Normalise v in situ - NB only shift non-zero digits */
	overflow = mpShiftLeft(v, v, shift, n);

	/* Copy normalised dividend u*d into r */
	overflow = mpShiftLeft(r, u, shift, n + m);
	uu = r;	/* Use ptr to keep notation constant */

	t[0] = overflow;	/* New digit Um+n */

	/* Step D2. Initialise j. Set j = m */
	for (j = m; j >= 0; j--)
	{
		/* Step D3. Calculate Qhat = (b.Uj+n + Uj+n-1)/Vn-1 */
		qhatOK = 0;
		t[1] = t[0];	/* This is Uj+n */
		t[0] = uu[j+n-1];
		overflow = spDivide(&qhat, &rhat, t, v[n-1]);

		/* Test Qhat */
		if (overflow)
		{	/* Qhat = b */
			qhat = VM_MAX_DIGIT;
			rhat = uu[j+n-1];
			rhat += v[n-1];
			if (rhat < v[n-1])	/* Overflow */
				qhatOK = 1;
		}
		if (!qhatOK && QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
		{	/* Qhat.Vn-2 > b.Rhat + Uj+n-2 */
			qhat--;
			rhat += v[n-1];
			if (!(rhat < v[n-1]))
				if (QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
					qhat--;
		}


		/* Step D4. Multiply and subtract */
		ww = &uu[j];
		overflow = mpMultSub(t[1], ww, v, qhat, (unsigned int)n);

		/* Step D5. Test remainder. Set Qj = Qhat */
		q[j] = qhat;
		if (overflow)
		{	/* Step D6. Add back if D4 was negative */
			q[j]--;
			overflow = mpAdd(ww, ww, v, (unsigned int)n);
		}

		t[0] = uu[j+n-1];	/* Uj+n on next round */

	}	/* Step D7. Loop on j */

	/* Clear high digits in uu */
	for (j = n; j < m+n; j++)
		uu[j] = 0;

	/* Step D8. Unnormalise. */

	mpShiftRight(r, r, shift, n);
	mpShiftRight(v, v, shift, n);

	return 0;
}

VM_DIGIT_T VMBigIntegers::mpMultSub(VM_DIGIT_T wn, VM_DIGIT_T w[], VM_DIGIT_T v[],
					   VM_DIGIT_T q, unsigned int n)
{	/*	Compute w = w - qv
		where w = (WnW[n-1]...W[0])
		return modified Wn.
	*/