mpilib.cpp

通訊保密編碼library project code.完整library project code!

/* C source code for multiprecision arithmetic library routines.
   Implemented Nov 86 by Philip Zimmermann
   Last revised 27 Nov 91 by PRZ

   Boulder Software Engineering
   3021 Eleventh Street
   Boulder, CO 80304
   [phone]

   (c) Copyright 1986-1994 by Philip Zimmermann.  All rights reserved.
   The author assumes no liability for damages resulting from the use
   of this software, even if the damage results from defects in this
   software.  No warranty is expressed or implied.  The use of this
   cryptographic software for developing weapon systems is expressly
   forbidden.


   These routines implement all of the multiprecision arithmetic
   necessary for number-theoretic cryptographic algorithms such as
   ElGamal, Diffie-Hellman, Rabin, or factoring studies for large
   composite numbers, as well as Rivest-Shamir-Adleman (RSA) public
   key cryptography.

   Although originally developed in Microsoft C for the IBM PC, this code
   contains few machine dependencies.  It assumes 2's complement
   arithmetic.  It can be adapted to 8-bit, 16-bit, or 32-bit machines,
   lowbyte-highbyte order or highbyte-lowbyte order.  This version
   has been converted to ANSI C.


   The internal representation for these extended precision integer
   "registers" is an array of "units".  A unit is a machine word, which
   is either an 8-bit byte, a 16-bit unsigned integer, or a 32-bit
   unsigned integer, depending on the machine's word size.  For example,
   an IBM PC or AT uses a unit size of 16 bits.  To perform arithmetic
   on these huge precision integers, we pass pointers to these unit
   arrays to various subroutines.  A pointer to an array of units is of
   type unitptr.  This is a pointer to a huge integer "register".

   When calling a subroutine, we always pass a pointer to the BEGINNING
   of the array of units, regardless of the byte order of the machine.
   On a lowbyte-first machine, such as the Intel 80x86, this unitptr
   points to the LEAST significant unit, and the array of units increases
   significance to the right.  On a highbyte-first machine, such as the
   Motorola 680x0, this unitptr points to the MOST significant unit, and
   the array of units decreases significance to the right.

   Modified 8 Apr 92 - HAJK
   Implement new VAX/VMS primitive support.

   Modified 30 Sep 92 -Castor Fu <[e-mail]>
   Upgraded PORTABLE support to allow sizeof(unit) == sizeof(long)

   Modified 28 Nov 92 - Thad Smith
   Added Smith modmult, generalized non-portable support.
 */

/* #define COUNTMULTS *//* count modmults for performance studies */

#ifdef DEBUG
#ifdef MSDOS
#ifdef __GO32__			/* DJGPP */
#include <pc.h>
#else
#include <conio.h>
#endif				/* __GO32__ */
#define poll_for_break() {while (kbhit()) getch();}
#endif				/* MSDOS */
#endif				/* DEBUG */

#ifndef poll_for_break
#define poll_for_break()	/* stub */
#endif

#include "mpilib.h"

#define min(a,b) (((a)<(b)) ? (a) : (b) )

#ifdef mp_smula
#ifdef mp_smul
Error:Both mp_smula and mp_smul cannot be defined.
#else
#define mp_smul	mp_smula
#endif
#endif

/* set macros for MULTUNIT */
#ifdef MUNIT8
#define MULTUNITSIZE   8
typedef unsigned char MULTUNIT;
#ifdef UNIT8
#define MULTUNIT_SIZE_SAME
#endif
#else				/* not MUNIT8 */
#ifdef MUNIT32
#define MULTUNITSIZE   32
typedef unsigned long MULTUNIT;
#ifdef UNIT32
#define MULTUNIT_SIZE_SAME
#else
/* #error is not portable, this has the same effect */
#include "UNITSIZE cannot be smaller than MULTUNITSIZE"
#endif
#else				/* assume MUNIT16 */
#define MULTUNITSIZE   16
typedef unsigned short MULTUNIT;
#ifdef UNIT16
#define MULTUNIT_SIZE_SAME
#endif				/* UNIT16 */
#ifdef UNIT8
#include "UNITSIZE cannot be smaller than MULTUNITSIZE"
#endif				/* UNIT8 */
#endif				/* MUNIT32 */
#endif				/* MUNIT8 */

#define MULTUNIT_msb    ((MULTUNIT)1 << (MULTUNITSIZE-1))	/* msb of
								   MULTUNIT */
#define DMULTUNIT_msb   (1L << (2*MULTUNITSIZE-1))
#define MULTUNIT_mask   ((MULTUNIT)((1L << MULTUNITSIZE)-1))
#define MULTUNITs_perunit   (UNITSIZE/MULTUNITSIZE)


void mp_smul(MULTUNIT * prod, MULTUNIT * multiplicand, MULTUNIT multiplier);
void mp_dmul(unitptr prod, const unit * multiplicand, const unit * multiplier);

short global_precision = MAX_UNIT_PRECISION;	/* units of precision for all routines */
/*      global_precision is the unit precision last set by set_precision.
   Initially, set_precision() should be called to define global_precision
   before using any of these other multiprecision library routines.
   i.e.:   set_precision(MAX_UNIT_PRECISION);
 */

/*************** multiprecision library primitives ****************/
/*      The following portable C primitives should be recoded in assembly.
   The entry point name should be defined, in "mpilib.h" to the external
   entry point name.  If undefined, the C version will be used.
 */

typedef unsigned long int ulint;

#ifndef mp_addc
/* 
   multiprecision add with carry r2 to r1, result in r1
   carry is incoming carry flag-- value should be 0 or 1
*/
boolean mp_addc
 (register unitptr r1, const unit * r2, register boolean carry)
{
    register unit x;
    short precision;		/* number of units to add */
    precision = global_precision;
    make_lsbptr(r1, precision);
    make_lsbptr(r2, precision);
    while (precision--) {
	if (carry) {
	    x = *r1 + *r2 + 1;
	    carry = (*r2 >= (unit) (~*r1));
	} else {
	    x = *r1 + *r2;
	    carry = (x < *r1);
	}
	post_higherunit(r2);
	*post_higherunit(r1) = x;
    }
    return carry;		/* return the final carry flag bit */
}				/* mp_addc */
#endif				/* mp_addc */

#ifndef mp_subb

/*
   multiprecision subtract with borrow, r2 from r1, result in r1
   borrow is incoming borrow flag-- value should be 0 or 1
 */
boolean mp_subb
 (register unitptr r1, const unit * r2, register boolean borrow)
{
    register unit x;
    short precision;		/* number of units to subtract */
    precision = global_precision;
    make_lsbptr(r1, precision);
    make_lsbptr(r2, precision);
    while (precision--) {
	if (borrow) {
	    x = *r1 - *r2 - 1;
	    borrow = (*r1 <= *r2);
	} else {
	    x = *r1 - *r2;
	    borrow = (*r1 < *r2);
	}
	post_higherunit(r2);
	*post_higherunit(r1) = x;
    }
    return borrow;		/* return the final carry/borrow flag bit */
}				/* mp_subb */
#endif	  			/* mp_subb */

#ifndef mp_rotate_left

/*
   multiprecision rotate left 1 bit with carry, result in r1.
   carry is incoming carry flag-- value should be 0 or 1
*/
boolean mp_rotate_left(register unitptr r1, register boolean carry)
{
    register int precision;	/* number of units to rotate */
    unsigned int mcarry = carry;/* int is supposed to be */
    unsigned int nextcarry=0;	/* the efficient size for ops */
	  					
    precision = global_precision;
    make_lsbptr(r1, precision);
    while (precision--) {
	nextcarry = (((signedunit) * r1) < 0);
	*r1 = (*r1 << 1) | mcarry;
	mcarry = nextcarry;
	pre_higherunit(r1);
    }
    return nextcarry;		/* return the final carry flag bit */
}				/* mp_rotate_left */
#endif				/* mp_rotate_left */

/************** end of primitives that should be in assembly *************/

/* The mp_shift_right_bits function is not called in any time-critical
   situations in public-key cryptographic functions, so it doesn't
   need to be coded in assembly language.
 */

/*
   multiprecision shift right bits, result in r1.
   bits is how many bits to shift, must be <= UNITSIZE.
*/
void mp_shift_right_bits(register unitptr r1, register short bits)
{
    register short precision;	/* number of units to shift */
    register unit carry, nextcarry, bitmask;
    register short unbits;
    if (bits == 0)
	return;			/* shift zero bits is a no-op */
    carry = 0;
    bitmask = power_of_2(bits) - 1;
    unbits = UNITSIZE - bits;	/* shift bits must be <= UNITSIZE */
    precision = global_precision;
    make_msbptr(r1, precision);
    if (bits == UNITSIZE) {
	while (precision--) {
	    nextcarry = *r1;
	    *r1 = carry;
	    carry = nextcarry;
	    pre_lowerunit(r1);
	}
    } else {
	while (precision--) {
	    nextcarry = *r1 & bitmask;
	    *r1 >>= bits;
	    *r1 |= carry << unbits;
	    carry = nextcarry;
	    pre_lowerunit(r1);
	}
    }
}				/* mp_shift_right_bits */

#ifndef mp_compare

/*
 * Compares multiprecision integers *r1, *r2, and returns:
 * -1 iff *r1 < *r2
 *  0 iff *r1 == *r2
 * +1 iff *r1 > *r2
 */
short mp_compare(const unit * r1, const unit * r2)
{
    register short precision;	/* number of units to compare */

    precision = global_precision;
    make_msbptr(r1, precision);
    make_msbptr(r2, precision);
    do {
	if (*r1 < *r2)
	    return -1;
	if (*post_lowerunit(r1) > *post_lowerunit(r2))
	    return 1;
    } while (--precision);
    return 0;			/*  *r1 == *r2  */
}				/* mp_compare */
#endif				/* mp_compare */

/* Increment multiprecision integer r. */
boolean mp_inc(register unitptr r)
{
    register short precision;
    precision = global_precision;
    make_lsbptr(r, precision);
    do {
	if (++(*r))
	    return FALSE;		/* no carry */
	post_higherunit(r);
    } while (--precision);
    return TRUE;			/* return carry set */
}				/* mp_inc */

/* Decrement multiprecision integer r. */
boolean mp_dec(register unitptr r)
{
    register short precision;
    precision = global_precision;
    make_lsbptr(r, precision);
    do {
	if ((signedunit) (--(*r)) != (signedunit) - 1)
	    return FALSE;		/* no borrow */
	post_higherunit(r);
    } while (--precision);
    return TRUE;			/* return borrow set */
}				/* mp_dec */

/* Compute 2's complement, the arithmetic negative, of r */
void mp_neg(register unitptr r)
{
    register short precision;	/* number of units to negate */
    precision = global_precision;
    mp_dec(r);			/* 2's complement is 1's complement plus 1 */
    do {			/* now do 1's complement */
	*r = ~(*r);
	r++;
    } while (--precision);
}				/* mp_neg */

#ifndef mp_move

void mp_move(register unitptr dst, register unitptr src)
{
    register short precision;	/* number of units to move */
    precision = global_precision;
    do {
	*dst++ = *src++;
    } while (--precision);
}				/* mp_move */
#endif				/* mp_move */

/* Init multiprecision register r with short value. */
void mp_init(register unitptr r, word16 value)
{	/* Note that mp_init doesn't extend sign bit for >32767 */

    unitfill0(r, global_precision);
    make_lsbptr(r, global_precision);
    *post_higherunit(r) = value;
#ifdef UNIT8
    *post_higherunit(r) = value >> UNITSIZE;
#endif
}				/* mp_init */

/* Returns number of significant units in r */
short significance(const unit * r)
{
    register short precision;
    precision = global_precision;
    make_msbptr(r, precision);
    do {
	if (*post_lowerunit(r))
	    return precision;
    } while (--precision);
    return precision;
}				/* significance */


#ifndef unitfill0

/* Zero-fill the unit buffer r. */
void unitfill0(unitptr r, word16 unitcount)
{
    while (unitcount--)
	*r++ = 0;
}				/* unitfill0 */
#endif				/* unitfill0 */

/* Unsigned divide, treats both operands as positive. */
int mp_udiv(register unitptr remainder, register unitptr quotient,
	    const unit * dividend, const unit * divisor)
{
    int bits;
    short dprec;
    register unit bitmask;

    if (testeq(divisor, 0))
	return -1;		/* zero divisor means divide error */
    mp_init0(remainder);