mpilib.cpp

各种加密算法的集合

/* 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 &amt; 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); 
mp_init0(quotient); 
/* normalize and compute number of bits in dividend first */ 
init_bitsniffer(dividend, bitmask, dprec, bits); 
/* rescale quotient to same precision (dprec) as dividend */ 
rescale(quotient, global_precision, dprec); 
make_msbptr(quotient, dprec); 

while (bits--) { 
mp_rotate_left(remainder, 
(boolean) (sniff_bit(dividend, bitmask) != 0));