vmbigintegersclass.cpp

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

/*****************************************************************************/
/*                              SOURCE FILE                                  */
/*****************************************************************************/
/*
       $Archive:   $

      $Revision:   $
          $Date:   $
        $Author:   $

    Description:   Implementation of big digits related functions. This 
                   collection of functions was packaged into this C++ class
                   in order to reduce the number of files needed to include
                   the functionality of the big numbers library. The original
                   authors 'header-comments' are shown below:

                   This source code is part of the BigDigits multiple-precision
                   arithmetic library Version 1.0 originally written by David Ireland,
                   copyright (c) 2001 D.I. Management Services Pty Limited, all rights
                   reserved. It is provided "as is" with no warranties. You may use
                   this software under the terms of the full copyright notice
                   "bigdigitsCopyright.txt" that should have been included with
                   this library. To obtain a copy send an email to
                   <[e-mail]> or visit <www.di-mgt.com.au/crypto.html>.
                   This notice must be retained in any copy.

                   See the header file for more copyright information.
*/
static char OBJECT_ID[] = "$Revision: 1 $ : $Date: 1/27/99 9:23a $";
/*****************************************************************************/

#include "../../../stdafx.h"
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <stdio.h>
#include "VMBigIntegersClass.h"


VM_DIGIT_T VM_SMALL_PRIMES_1[] = 
{ 
  2, 3, 5, 7, 11, 13, 17, 19 
};
#define VM_N_SMALL_PRIMES_1 sizeof( VM_SMALL_PRIMES_1 ) / sizeof( VM_DIGIT_T )


#define B ( VM_MAX_HALF_DIGIT + 1 )

VM_DIGIT_T VM_SMALL_PRIMES_2[] = 
{ 
   3,    5,   7,   11,   13,   17, 
  19,   23,  29,   31,   37,   41, 
  43,   47,  53,   59,   61,   67, 
  71,   73,  79,   83,   89,   97, 
  101, 103, 107 
};
#define VM_N_SMALL_PRIMES_2 sizeof(VM_SMALL_PRIMES_2)/sizeof(VM_DIGIT_T)


//char *vm_copyright_notice(void)
//{
//	return 
//"Contains multiple-precision arithmetic code originally written by David Ireland, copyright (c) 2001 by D.I. Management Services Pty Limited <www.di-mgt.com.au>, and is used with permission.";
//}



VM_DIGIT_T VMBigIntegers::spPseudoRand( VM_DIGIT_T lower, VM_DIGIT_T upper )
{	/*	Returns a pseudo-random digit.
		Handles own seeding using time
		NOT for cryptographically-secure random numbers.
	*/
	static unsigned seeded = 0;
	unsigned int i;
	double f;

	if (!seeded)
	{
		/* seed with system time */
		srand((unsigned)time(NULL));
		/* Throw away a random few to avoid similar starts */
		i = rand() & 0xFF;
		while (i--)
			rand();
		seeded = 1;
	}
	f = (double)rand() / RAND_MAX * upper;
	return (lower + (VM_DIGIT_T)f);
}


int VMBigIntegers::spMultiply(VM_DIGIT_T p[2], VM_DIGIT_T x, VM_DIGIT_T y)
{	/*	Computes p = x * y */
	/*	Ref: Arbitrary Precision Computation
	http://numbers.computation.free.fr/Constants/constants.html

         high    p1                p0     low
        +--------+--------+--------+--------+
        |      x1*y1      |      x0*y0      |
        +--------+--------+--------+--------+
               +-+--------+--------+
               |1| (x0*y1 + x1*y1) |
               +-+--------+--------+
                ^carry from adding (x0*y1+x1*y1) together
                        +-+
                        |1|< carry from adding VM_LOHALF t
                        +-+  to high half of p0
	*/
	VM_DIGIT_T x0, y0, x1, y1;
	VM_DIGIT_T t, u, carry;

	/*	Split each x,y into two halves
		x = x0 + B*x1
		y = y0 + B*y1
		where B = 2^16, half the digit size
		Product is
		xy = x0y0 + B(x0y1 + x1y0) + B^2(x1y1)
	*/

	x0 = VM_LOHALF(x);
	x1 = VM_HIHALF(x);
	y0 = VM_LOHALF(y);
	y1 = VM_HIHALF(y);

	/* Calc low part - no carry */
	p[0] = x0 * y0;

	/* Calc middle part */
	t = x0 * y1;
	u = x1 * y0;
	t += u;
	if (t < u)
		carry = 1;
	else
		carry = 0;

	/*	This carry will go to high half of p[1]
		+ high half of t into low half of p[1] */
	carry = VM_TOHIGH(carry) + VM_HIHALF(t);

	/* Add low half of t to high half of p[0] */
	t = VM_TOHIGH(t);
	p[0] += t;
	if (p[0] < t)
		carry++;

	p[1] = x1 * y1;
	p[1] += carry;


	return 0;
}


int VMBigIntegers::spModMult(VM_DIGIT_T *a, VM_DIGIT_T x, VM_DIGIT_T y, VM_DIGIT_T m)
{	/*	Computes a = (x * y) mod m */
	
	/* Calc p[2] = x * y */
	VM_DIGIT_T p[2];
	spMultiply(p, x, y);

	/* Then modulo */
	*a = mpShortMod(p, m, 2);
	return 0;
}


int VMBigIntegers::spModInv(VM_DIGIT_T *inv, VM_DIGIT_T u, VM_DIGIT_T v)
{	/*	Computes inv = u^(-1) mod v */
	/*	Ref: Knuth Algorithm X Vol 2 p 342 
		ignoring u2, v2, t2
		and avoiding negative numbers
	*/
	VM_DIGIT_T u1, u3, v1, v3, t1, t3, q, w;
	int bIterations = 1;
	
	/* Step X1. Initialise */
	u1 = 1;
	u3 = u;
	v1 = 0;
	v3 = v;

	while (v3 != 0)	/* Step X2. */
	{	/* Step X3. */
		q = u3 / v3;	/* Divide and */
		t3 = u3 % v3;
		w = q * v1;		/* "Subtract" */
		t1 = u1 + w;
		/* Swap */
		u1 = v1;
		v1 = t1;
		u3 = v3;
		v3 = t3;
		bIterations = -bIterations;
	}

	if (bIterations < 0)
		*inv = v - u1;
	else
		*inv = u1;

	return 0;
}



int VMBigIntegers::spModExp(VM_DIGIT_T *exp, VM_DIGIT_T x,  VM_DIGIT_T n, VM_DIGIT_T d )
{
	return spModExpB(exp, x, n, d);
}


int VMBigIntegers::spModExpK(VM_DIGIT_T *exp, VM_DIGIT_T x, VM_DIGIT_T n, VM_DIGIT_T d )
{	/*	Computes exp = x^n mod d */
	/*	Ref: Knuth Vol 2 Ch 4.6.3 p 462 Algorithm A
	*/
	VM_DIGIT_T y = 1;		/* Step A1. Initialise */

	while (n > 0)
	{							/* Step A2. Halve N */
		if (n & 0x1)			/* If odd */
			spModMult(&y, y, x, d);	/*   Step A3. Multiply Y by Z */	
		
		n >>= 1;					/* Halve N */
		if (n > 0)				/* Step A4. N = 0? Y is answer */
			spModMult(&x, x, x, d);	/* Step A5. Square Z */
	}

	*exp = y;
	return 0;
}

int VMBigIntegers::spModExpB(VM_DIGIT_T *exp, VM_DIGIT_T x, 
			VM_DIGIT_T e, VM_DIGIT_T m)
{	/*	Computes exp = x^e mod m */
	/*	Binary left-to-right method
	*/
	VM_DIGIT_T mask;
	VM_DIGIT_T y;	/* Temp variable */

	/* Find most significant bit in e */
	for (mask = VM_HIBITMASK; mask > 0; mask >>= 1)
	{
		if (e & mask)
			break;
	}

	y = x;
	/* For j = k-2 downto 0 step -1 */
	for (mask >>= 1; mask > 0; mask >>= 1)
	{
		spModMult(&y, y, y, m);		/* Square */
		if (e & mask)
			spModMult(&y, y, x, m);	/* Multiply */
	}

	*exp = y;
	return 0;
}


int VMBigIntegers::spIsPrime(VM_DIGIT_T w, 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.
	*/
	/*	Rabin-Miller Probabilistic Primality Test,
		from FIPS-186-2 Appendix 2.
		Also Schneier 2nd ed p 260 & Knuth Vol 2, p 379.
	*/

	unsigned int i, j;
	VM_DIGIT_T m, a, b, z;
	int failed;

	/*	First check for small primes */
	for (i = 0; i < VM_N_SMALL_PRIMES_1; i++)
	{
		if (w % VM_SMALL_PRIMES_1[i] == 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 */
	m = w - 1;
	for (a = 0; VM_ISEVEN(m); a++)
		m >>= 1;	/* Divide by 2 until m is odd */

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

	for (i = 0; i < t; i++)
	{
		failed = 1;	/* Assume fail unless passed in loop */
		/* Step 3. Generate a random integer 1 < b < w */
		b = spPseudoRand(2, w - 1);

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

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

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

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

		if (failed)
		{	/* Step 8. Not a prime - stop */
			return 0;
		}
	}	/* Step 9. Go to Step 3 until i >= n */
	/* If got here, probably prime => success */
	return 1;
}


VM_DIGIT_T VMBigIntegers::spGcd(VM_DIGIT_T x, VM_DIGIT_T y)
{	/*	Returns gcd(x, y) */

	/* Ref: Schneier 2nd ed, p245 */
	
	VM_DIGIT_T g;

	if (x + y == 0)
		return 0;	/* Error */

	g = y;
	while (x > 0)
	{
		g = x;
		x = y % x;
		y = g;
	}
	return g;
}


VM_DIGIT_T VMBigIntegers::spDivide(VM_DIGIT_T *q, VM_DIGIT_T *r, VM_DIGIT_T u[2], VM_DIGIT_T v)
{	/*	Computes quotient q = u / v, remainder r = u mod v
		where u is a double digit
		and q, v, r are single precision digits.
		Returns high digit of quotient (max value is 1)
		Assumes normalised such that v1 >= b/2
		where b is size of HALF_DIGIT
		i.e. the most significant bit of v should be one

		In terms of half-digits in Knuth notation:
		(q2q1q0) = (u4u3u2u1u0) / (v1v0)
		(r1r0) = (u4u3u2u1u0) mod (v1v0)
		for m = 2, n = 2 where u4 = 0
		q2 is either 0 or 1.
		We set q = (q1q0) and return q2 as "overflow'
	*/
	VM_DIGIT_T qhat, rhat, t, v0, v1, u0, u1, u2, u3;
	VM_DIGIT_T uu[2], q2;

	/* Check for normalisation */
	if (!(v & VM_HIBITMASK))
	{
		*q = *r = 0;
		return VM_MAX_DIGIT;
	}
	
	/* Split up into half-digits */
	v0 = VM_LOHALF(v);
	v1 = VM_HIHALF(v);
	u0 = VM_LOHALF(u[0]);
	u1 = VM_HIHALF(u[0]);
	u2 = VM_LOHALF(u[1]);
	u3 = VM_HIHALF(u[1]);

	/* Do three rounds of Knuth Algorithm D Vol 2 p272 */

	/*	ROUND 1. Set j = 2 and calculate q2 */
	/*	Estimate qhat = (u4u3)/v1  = 0 or 1 
		then set (u4u3u2) -= qhat(v1v0)
		where u4 = 0.
	*/
	qhat = u3 / v1;
	if (qhat > 0)
	{
		rhat = u3 - qhat * v1;
		t = VM_TOHIGH(rhat) | u2;
		if (qhat * v0 > t)
			qhat--;
	}
	uu[1] = 0;		/* (u4) */
	uu[0] = u[1];	/* (u3u2) */
	if (qhat > 0)
	{
		/* (u4u3u2) -= qhat(v1v0) where u4 = 0 */
		spMultSub(uu, qhat, v1, v0);
		if (VM_HIHALF(uu[1]) != 0)
		{	/* Add back */
			qhat--;
			uu[0] += v;
			uu[1] = 0;
		}
	}
	q2 = qhat;

	/*	ROUND 2. Set j = 1 and calculate q1 */
	/*	Estimate qhat = (u3u2) / v1 
		then set (u3u2u1) -= qhat(v1v0)
	*/
	t = uu[0];
	qhat = t / v1;
	rhat = t - qhat * v1;
	/* Test on v0 */
	t = VM_TOHIGH(rhat) | u1;
	if ((qhat == B) || (qhat * v0 > t))
	{
		qhat--;
		rhat += v1;
		t = VM_TOHIGH(rhat) | u1;
		if ((rhat < B) && (qhat * v0 > t))
			qhat--;
	}

	/*	Multiply and subtract 
		(u3u2u1)' = (u3u2u1) - qhat(v1v0)	
	*/
	uu[1] = VM_HIHALF(uu[0]);	/* (0u3) */
	uu[0] = VM_TOHIGH(VM_LOHALF(uu[0])) | u1;	/* (u2u1) */
	spMultSub(uu, qhat, v1, v0);
	if (VM_HIHALF(uu[1]) != 0)
	{	/* Add back */
		qhat--;
		uu[0] += v;
		uu[1] = 0;
	}

	/* q1 = qhat */
	*q = VM_TOHIGH(qhat);

	/* ROUND 3. Set j = 0 and calculate q0 */
	/*	Estimate qhat = (u2u1) / v1
		then set (u2u1u0) -= qhat(v1v0)
	*/
	t = uu[0];
	qhat = t / v1;
	rhat = t - qhat * v1;
	/* Test on v0 */
	t = VM_TOHIGH(rhat) | u0;
	if ((qhat == B) || (qhat * v0 > t))
	{
		qhat--;
		rhat += v1;
		t = VM_TOHIGH(rhat) | u0;
		if ((rhat < B) && (qhat * v0 > t))
			qhat--;
	}

	/*	Multiply and subtract 
		(u2u1u0)" = (u2u1u0)' - qhat(v1v0)
	*/
	uu[1] = VM_HIHALF(uu[0]);	/* (0u2) */
	uu[0] = VM_TOHIGH(VM_LOHALF(uu[0])) | u0;	/* (u1u0) */
	spMultSub(uu, qhat, v1, v0);
	if (VM_HIHALF(uu[1]) != 0)
	{	/* Add back */
		qhat--;
		uu[0] += v;
		uu[1] = 0;
	}

	/* q0 = qhat */
	*q |= VM_LOHALF(qhat);

	/* Remainder is in (u1u0) i.e. uu[0] */
	*r = uu[0];
	return q2;
}


void VMBigIntegers::spMultSub(VM_DIGIT_T uu[2], VM_DIGIT_T qhat, 
					  VM_DIGIT_T v1, VM_DIGIT_T v0)
{
	/*	Compute uu = uu - q(v1v0) 
		where uu = u3u2u1u0, u3 = 0
		and u_n, v_n are all half-digits
		even though v1, v2 are passed as full digits.
	*/
	VM_DIGIT_T p0, p1, t;

	p0 = qhat * v0;