t_mptest.c

大数的计算包括加减乘除

/* t_mpTest.c */

/******************* SHORT COPYRIGHT NOTICE*************************
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.
****************** END OF COPYRIGHT NOTICE*************************/

#include <stdio.h>
#include <assert.h>
#include "bigdigits.h"

#define TEST_LEN MAX_DIG_LEN	

DIGIT_T w[TEST_LEN]; 
DIGIT_T u[TEST_LEN];
DIGIT_T v[TEST_LEN];
DIGIT_T q[TEST_LEN];
DIGIT_T r[TEST_LEN];
DIGIT_T p1[TEST_LEN];
DIGIT_T p2[TEST_LEN];
DIGIT_T p3[TEST_LEN];
DIGIT_T g[TEST_LEN];

void ClearAll(void)
{
	mpSetZero(u, TEST_LEN);
	mpSetZero(v, TEST_LEN);
	mpSetZero(w, TEST_LEN);
	mpSetZero(q, TEST_LEN);
	mpSetZero(r, TEST_LEN);
	mpSetZero(g, TEST_LEN);
	mpSetZero(p1, TEST_LEN);
	mpSetZero(p2, TEST_LEN);
	mpSetZero(p3, TEST_LEN);
}

static int MakeMultiplePrime(DIGIT_T p[], unsigned int ndigits)
{	/* Returns a random prime number of ndigits */
	/*	WARNING: This is not cryptographically secure 
		because the random number generator isn't */
	unsigned int i;

	for (i = 0; i < ndigits; i++)
		p[i] = spPseudoRand(0, MAX_DIGIT);

	/* Make sure the highest and low bits are set */
	p[ndigits - 1] |= HIBITMASK;
	p[0] |= 0x1;

	/*printf("p="); mpPrintNL(p, ndigits);*/

	/* Check if prime */
	while (!mpIsPrime(p, ndigits, 10))
	{
		/* Keep adding 2 until find a prime */
		mpShortAdd(p, p, 2, ndigits);

		/*printf("p="); mpPrintNL(p, ndigits);*/
		printf(".");

		/* Check for overflow */
		if (!(p[ndigits - 1] & HIBITMASK))
			return -1;	/* Failed to find a prime */
	}

	return 0;
}

unsigned int MakeMultipleRandom(DIGIT_T a[], unsigned int ndigits)
{
	/* Make a random number of up to ndigits digits */ 
	unsigned int i, n, bits;
	DIGIT_T mask;

	n = (unsigned int)spPseudoRand(1, ndigits);
	for (i = 0; i < n; i++)
		a[i] = spPseudoRand(0, MAX_DIGIT);
	for (i = n; i < ndigits; i++)
		a[i] = 0;

	/*	Zero out a random number of bits in leading digit 
		about half the time */
	bits = (unsigned int)spPseudoRand(0, 2*BITS_PER_DIGIT);
	if (bits != 0 && bits < BITS_PER_DIGIT)
	{
		mask = HIBITMASK;
		for (i = 1; i < bits; i++)
		{
			mask |= (mask >> 1);
		}
		mask = ~mask;
		a[n-1] &= mask;
	}
	return n;
}

void ShowAdd(DIGIT_T w[], DIGIT_T u[], DIGIT_T v[], 
			 DIGIT_T carry, unsigned int ndigits)
{
	printf("mpAdd: ");
	mpPrintTrim(u, ndigits);
	printf("+ ");
	mpPrintTrim(v, ndigits);
	printf("= ");
	mpPrintTrim(w, ndigits);
	printf(", Carry = %lx\n", carry);
}

main()
{
	DIGIT_T carry, m;

	/*	Force linker to include copyright notice in 
		executable object image
	*/
	copyright_notice();
		
	ClearAll();

	/* Start easy: 1 + 1 = 2 */
	mpSetDigit(u, 1, TEST_LEN);	/* u = 1 */
	carry = mpAdd(w, u, u, TEST_LEN);	/* w = u + u */
	ShowAdd(w, u, u, carry, TEST_LEN);
	/* Check that w == 2 */
	assert(mpShortCmp(w, 2, TEST_LEN) == 0);

	/* ---- Add and subtract ---- */
	/* Add two random numbers w = u + v */
	MakeMultipleRandom(u, TEST_LEN-1);
	MakeMultipleRandom(v, TEST_LEN-1);
	carry = mpAdd(w, u, v, TEST_LEN);

	/* r = w - v */
	carry = mpSubtract(r, w, v, TEST_LEN);

	/* Check r == u */
	assert(mpEqual(r, u, TEST_LEN));
	printf("Add and subtract worked OK\n");
	ClearAll();

	/* ---- Multiply and divide ---- */
	/* Multiply two random numbers w = u * v */
	MakeMultipleRandom(u, TEST_LEN / 2);
	MakeMultipleRandom(v, TEST_LEN / 2);
	mpMultiply(w, u, v, TEST_LEN / 2);

	/* q = w / v, r = w % v */
	mpDivide(q, r, w, TEST_LEN, v, TEST_LEN / 2);
	/* Check q == u and r == 0 */
	assert(mpEqual(q, u, TEST_LEN / 2));
	assert(mpIsZero(r, TEST_LEN / 2));

	ClearAll();
	
	/* Pick two random numbers u, v */
	MakeMultipleRandom(u, TEST_LEN);
	MakeMultipleRandom(v, TEST_LEN);
	/* Divide one by the other: q = u / v, r = u % v */
	mpDivide(q, r, u, TEST_LEN, v, TEST_LEN);
	/* Check w = q * v + r == u */
	mpMultiply(w, q, v, TEST_LEN);
	mpAdd(w, w, r, TEST_LEN);
	assert(mpEqual(w, u, TEST_LEN));
	printf("Multiply and divide worked OK\n");

	/* ---- Greatest Common Divisor ---- */
	/* Pick 3x random primes p1, p2, p3 */
	printf("Creating 3 prime numbers (be patient):\n");
	MakeMultiplePrime(p1, TEST_LEN / 2);
	printf("(1)\n");
	MakeMultiplePrime(p2, TEST_LEN / 2);
	printf("(2)\n");
	MakeMultiplePrime(p3, TEST_LEN / 2);
	printf("(3)\n");

	/* Calculate two products from these primes */
	/* u = p1 * p2 */
	mpMultiply(u, p1, p2, TEST_LEN / 2);
	/* v = p1 * p3 */
	mpMultiply(v, p1, p3, TEST_LEN / 2);
	/* Now calculate g = gcd(u, v) */
	mpGcd(g, u, v, TEST_LEN);
	/* And check that g == p1 */
	assert(mpEqual(g, p1, TEST_LEN)); 
	printf("Greatest Common Divisor worked OK\n");

	/* ---- Modular Inverse ---- */
	/* Use previous prime as modulus, v */
	mpSetEqual(v, p1, TEST_LEN);

	/* Pick a small multiplier, m */
	m = spPseudoRand(1, MAX_DIGIT);

	/* Set u = (vm - 1) */
	mpShortMult(u, v, m, TEST_LEN);
	mpShortSub(u, u, 1, TEST_LEN);

	mpModInv(w, u, v, TEST_LEN);
	/* Check that g = (w * u) mod v = 1 */
	mpModMult(g, w, u, v, TEST_LEN);
	assert((mpShortCmp(g, 1, TEST_LEN) == 0));
	printf("Modular inverse worked OK\n");

	/* For further checks do RSA calculation - see t_mpRSA.c */
	
	return 0;
}