mpilib.cpp

各种加密算法的集合

*/ 
mp = modulus; 
sigmod = significance(modulus); 
rescale(mp, original_precision, sigmod); 
/* prec is the unit precision required for 3 MULTUNITs */ 
prec = (3 + (MULTUNITs_perunit - 1)) / MULTUNITs_perunit; 
set_precision(prec); 

/* set mp = ptr to most significant units of modulus, then move 
* the most significant units to mp_divisor 
*/ 
mp = msbptr(mp, sigmod) - tohigher(prec - 1); 
unmake_lsbptr(mp, prec); 
mp_move(mod_divisor, mp); 

/* Keep 2*MULTUNITSIZE bits in mod_divisor. 
* This will (normally) result in a reciprocal of 2*MULTUNITSIZE+1 bits. 
*/ 
mshift = countbits(mod_divisor) - 2 * MULTUNITSIZE; 
mp_shift_right_bits(mod_divisor, mshift); 
mp_recip(mod_quotient, mod_divisor); 
mp_shift_right_bits(mod_quotient, 1); 

/* Reduce to: 0 < mshift <= MULTUNITSIZE */ 
mshift = ((mshift + (MULTUNITSIZE - 1)) > MULTUNITSIZE) + 1; 
/* round up, rescaling if necessary */ 
mp_inc(mod_quotient); 
if (countbits(mod_quotient) > 2 * MULTUNITSIZE) { 
mp_shift_right_bits(mod_quotient, 1); 
mshift--; /* now 0 <= mshift <= MULTUNITSIZE */ 
} 
mpm = lsbptr((MULTUNIT *) mod_quotient, prec * MULTUNITs_perunit); 
recipl = *post_higherunit(mpm); 
reciph = *mpm; 
mp_set_recip(reciph, recipl, mshift); 
set_precision(original_precision); 
return 0; /* normal return */ 
} /* stage_smith_modulus */ 

/* Smith's algorithm performs a multiply combined with a modulo operation. 
Computes: prod = (multiplicand*multiplier) mod modulus 
The modulus must first be set by stage_smith_modulus(). 
smith_modmult() is aliased to mp_modmult(). 
*/ 
int smith_modmult(unitptr prod, const unit * multiplicand, const unit * multiplier) 
{ 
unitptr d; /* ptr to product */ 
MULTUNIT *dmph, *dmpl, *dmp; /* ptrs to dividend (high, low, first) 
* aligned for subtraction */ 
/* Note that dmph points one MULTUNIT higher than indicated by 
global precision. This allows us to zero out a word one higher than 
the normal precision. 
*/ 
short orig_precision; 
short nqd; /* number of quotient digits remaining to 
* be generated */ 
short dmi; /* number of significant MULTUNITs in 
product */ 

d = ds_data + 1; /* room for leading MSB if HIGHFIRST */ 
orig_precision = global_precision; 
mp_dmul(d, multiplicand, multiplier); 

rescale(d, orig_precision * 2, orig_precision * 2 + 1); 
set_precision(orig_precision * 2 + 1); 
*msbptr(d, global_precision) = 0; /* leading 0 unit */ 

/* We now start working with MULTUNITs. 
Determine the most significant MULTUNIT of the product so we don't 
have to process leading zeros in our divide loop. 
*/ 
dmi = significance(d) * MULTUNITs_perunit; 
if (dmi >= modlenMULTUNITS) { 
/* Make dividend negative. This allows the use of mp_smul to 
* "subtract" the product of the modulus and the trial divisor 
* by actually adding to a negative dividend. 
* The one's complement of the dividend is used, since it causes 
* a zero value to be represented as all ones. This facilitates 
* testing the result for possible overflow, since a sign bit 
* indicates that no adjustment is necessary, and we should not 
* attempt to adjust if the result of the addition is zero. 
*/ 
mp_inc(d); 
mp_neg(d); 
set_precision(orig_precision); 
munit_prec = global_precision * UNITSIZE / MULTUNITSIZE; 

/* Set msb, lsb, and normal ptrs of dividend */ 
dmph = lsbptr((MULTUNIT *) d, (orig_precision * 2 + 1) * 
MULTUNITs_perunit) + tohigher(dmi); 
nqd = dmi + 1 - modlenMULTUNITS; 
dmpl = dmph - tohigher(modlenMULTUNITS); 

/* 
* Divide loop. 
* Each iteration computes the next quotient MULTUNIT digit, then 
* multiplies the divisor (modulus) by the quotient digit and adds 
* it to the one's complement of the dividend (equivalent to 
* subtracting). If the product was greater than the remaining dividend, 
* we get a non-negative result, in which case we subtract off the 
* modulus to get the proper negative remainder. 
*/ 
for (; nqd; nqd--) { 
MULTUNIT q; /* quotient trial digit */ 

q = mp_quo_digit(dmph); 
if (q > 0) { 
mp_smul(dmpl, modmpl, q); 

/* Perform correction if q too large. 
* This rarely occurs. 
*/ 
if (!(*dmph &amt; MULTUNIT_msb)) { 
dmp = dmpl; 
unmake_lsbptr(dmp, orig_precision * 
MULTUNITs_perunit); 
if (mp_musubb(dmp, 
(MULTUNIT *) modulus, FALSE)) 
(*dmph)--; 
} 
} 
pre_lowerunit(dmph); 
pre_lowerunit(dmpl); 
} 
/* d contains the one's complement of the remainder. */ 
rescale(d, orig_precision * 2 + 1, orig_precision); 
set_precision(orig_precision); 
mp_neg(d); 
mp_dec(d); 
} else { 
/* Product was less than modulus. Return it. */ 
rescale(d, orig_precision * 2 + 1, orig_precision); 
set_precision(orig_precision); 
} 
mp_move(prod, d); 
return (0); /* normal return */ 
} /* smith_modmult */ 

/* Smith's mp_modmult function leaves some internal arrays in memory, 
so we have to call modmult_burn() at the end of mp_modexp. 
This is so that no cryptographically sensitive data is left in memory 
after the program exits. 
smith_burn() is aliased to modmult_burn(). 
*/ 
void smith_burn(void) 
{ 
empty_array(modulus); 
empty_array(ds_data); 
empty_array(mod_quotient); 
empty_array(mod_divisor); 
modmpl = 0; 
mshift = nbits = 0; 
reciph = recipl = 0; 
modlenMULTUNITS = mutemp = 0; 
mp_set_recip(0, 0, 0); 
} /* smith_burn */ 

/* End of Thad Smith's implementation of modmult. */ 

#endif /* SMITH */ 

/* Returns number of significant bits in r */ 
int countbits(const unit * r) 
{ 
int bits; 
short prec; 
register unit bitmask; 
init_bitsniffer(r, bitmask, prec, bits); 
return bits; 
} /* countbits */ 

char *copyright_notice(void) 
/* force linker to include copyright notice in the executable object image. */ 
{ 
return ("(c)1986 Philip Zimmermann"); 
} /* copyright_notice */ 

/* 
* Russian peasant combined exponentiation/modulo algorithm. 
* Calls modmult instead of mult. 
* Computes: expout = (expin**exponent) mod modulus 
* WARNING: All the arguments must be less than the modulus! 
*/ 
int mp_modexp(register unitptr expout, const unit * expin, 
const unit * exponent, const unit * modulus) 
{ 
int bits; 
short oldprecision; 
register unit bitmask; 
unit product[MAX_UNIT_PRECISION]; 
short eprec; 

#ifdef COUNTMULTS 
tally_modmults = 0; /* clear "number of modmults" counter */ 
tally_modsquares = 0; /* clear "number of modsquares" counter */ 
#endif /* COUNTMULTS */ 
mp_init(expout, 1); 
if (testeq(exponent, 0)) { 
if (testeq(expin, 0)) 
return -1; /* 0 to the 0th power means return error */ 
return 0; /* otherwise, zero exponent means expout 
is 1 */ 
} 
if (testeq(modulus, 0)) 
return -2; /* zero modulus means error */ 
#if SLOP_BITS > 0 /* if there's room for sign bits */ 
if (mp_tstminus(modulus)) 
return -2; /* negative modulus means error */ 
#endif /* SLOP_BITS > 0 */ 
if (mp_compare(expin, modulus) >= 0) 
return -3; /* if expin >= modulus, return error */ 
if (mp_compare(exponent, modulus) >= 0) 
return -4; /* if exponent >= modulus, return error */ 

oldprecision = global_precision; /* save global_precision */ 
/* set smallest optimum precision for this modulus */ 
set_precision(bits2units(countbits(modulus) + SLOP_BITS)); 
/* rescale all these registers to global_precision we just defined */ 
rescale(modulus, oldprecision, global_precision); 
rescale(expin, oldprecision, global_precision); 
rescale(exponent, oldprecision, global_precision); 
rescale(expout, oldprecision, global_precision); 

if (stage_modulus(modulus)) { 
set_precision(oldprecision); /* restore original precision */ 
return -5; /* unstageable modulus (STEWART algorithm) */ 
} 
/* normalize and compute number of bits in exponent first */ 
init_bitsniffer(exponent, bitmask, eprec, bits); 

/* We can "optimize out" the first modsquare and modmult: */ 
bits--; /* We know for sure at this point that 
bits>0 */ 
mp_move(expout, expin); /* expout = (1*1)*expin; */ 
bump_bitsniffer(exponent, bitmask); 

while (bits--) { 
poll_for_break(); /* polls keyboard, allows ctrl-C to 
abort program */ 
#ifdef COUNTMULTS 
tally_modsquares++; /* bump "number of modsquares" counter */ 
#endif /* COUNTMULTS */ 
mp_modsquare(product, expout); 
if (sniff_bit(exponent, bitmask)) { 
mp_modmult(expout, product, expin); 
#ifdef COUNTMULTS 
tally_modmults++; /* bump "number of modmults" counter */ 
#endif /* COUNTMULTS */ 
} else { 
mp_move(expout, product); 
} 
bump_bitsniffer(exponent, bitmask); 
} /* while bits-- */ 
mp_burn(product); /* burn the evidence on the stack */ 
modmult_burn(); /* ask mp_modmult to also burn its 
own evidence */ 

#ifdef COUNTMULTS /* diagnostic analysis */ 
{ 
long atomic_mults; 
unsigned int unitcount, totalmults; 
unitcount = bits2units(countbits(modulus)); 
/* calculation assumes modsquare takes as long as a modmult: */ 
atomic_mults = (long) tally_modmults *(unitcount * unitcount); 
atomic_mults += (long) tally_modsquares *(unitcount * unitcount); 
printf(">ld atomic mults for ", atomic_mults); 
printf(">d+>d = >d modsqr+modmlt, at >d bits, >d words.\n", 
tally_modsquares, tally_modmults, 
tally_modsquares + tally_modmults, 
countbits(modulus), unitcount); 
} 
#endif /* COUNTMULTS */ 

set_precision(oldprecision); /* restore original precision */ 

/* Do an explicit reference to the copyright notice so that the linker 
will be forced to include it in the executable object image... */ 
copyright_notice(); /* has no real effect at run time */ 

return 0; /* normal return */ 
} /* mp_modexp */ 

/* 
* This is a faster modexp for moduli with a known factorisation into two 
* relatively prime factors p and q, and an input relatively prime to the 
* modulus, the Chinese Remainder Theorem to do the computation mod p and 
* mod q, and then combine the results. This relies on a number of 
* precomputed values, but does not actually require the modulus n or the 
* exponent e. 
* 
* expout = expin ^ e mod (p*q). 
* We form this by evaluating 
* p2 = (expin ^ e) mod p and 
* q2 = (expin ^ e) mod q 
* and then combining the two by the CRT. 
* 
* Two optimisations of this are possible. First, we can reduce expin 
* modulo p and q before starting. 
* 
* Second, since we know the factorisation of p and q (trivially derived 
* from the factorisation of n = p*q), and expin is relatively prime to 
* both p and q, we can use Euler's theorem, expin^phi(m) = 1 (mod m), 
* to throw away multiples of phi(p) or phi(q) in e. 
* Letting ep = e mod phi(p) and 
* eq = e mod phi(q) 
* then combining these two speedups, we only need to evaluate 
* p2 = ((expin mod p) ^ ep) mod p and 
* q2 = ((expin mod q) ^ eq) mod q. 
* 
* Now we need to apply the CRT. Starting with 
* expout = p2 (mod p) and 
* expout = q2 (mod q) 
* we can say that expout = p2 + p * k, and if we assume that 0 <= p2 < p, 
* then 0 <= expout < p*q for some 0 <= k < q. Since we want expout = q2 
* (mod q), then p*k = q2-p2 (mod q). Since p and q are relatively prime, 
* p has a multiplicative inverse u mod q. In other words, u = 1/p (mod q). 
* 
* Multiplying by u on both sides gives k = u*(q2-p2) (mod q). 
* Since we want 0 <= k < q, we can thus find k as 
* k = (u * (q2-p2)) mod q. 
* 
* Once we have k, evaluating p2 + p * k is easy, and 
* that gives us the result. 
* 
* In the detailed implementation, there is a temporary, temp, used to 
* hold intermediate results, p2 is held in expout, and q2 is used as a 
* temporary in the final step when it is no longer needed. With that, 
* you should be able to understand the code below. 
*/ 
int mp_modexp_crt(unitptr expout, const unit * expin, 
const unit * p, const unit * q, const unit * ep, const unit * eq, const unit * u) 
{ 
unit q2[MAX_UNIT_PRECISION]; 
unit temp[MAX_UNIT_PRECISION]; 
int status; 

/* First, compiute p2 (physically held in M) */ 

/* p2 = [ (expin mod p) ^ ep ] mod p */ 
mp_mod(temp, expin, p); /* temp = expin mod p */ 
status = mp_modexp(expout, temp, ep, p); 
if (status < 0) { /* mp_modexp returned an error. */ 
mp_init(expout, 1); 
return status; /* error return */ 
} 
/* And the same thing for q2 */ 

/* q2 = [ (expin mod q) ^ eq ] mod q */ 
mp_mod(temp, expin, q); /* temp = expin mod q */ 
status = mp_modexp(q2, temp, eq, q); 
if (status < 0) { /* mp_modexp returned an error. */ 
mp_init(expout, 1); 
return status; /* error return */ 
} 
/* Now use the multiplicative inverse u to glue together the 
two halves. 
*/ 
#if 0 
/* This optimisation is useful if you commonly get small results, 
but for random results and large q, the odds of (1/q) of it 
being useful do not warrant the test. 
*/ 
if (mp_compare(expout, q2) != 0) { 
#endif 
/* Find q2-p2 mod q */ 
if (mp_sub(q2, expout)) /* if the result went negative */ 
mp_add(q2, q); /* add q to q2 */ 

/* expout = p2 + ( p * [(q2*u) mod q] ) */ 
mp_mult(temp, q2, u); /* q2*u */ 
mp_mod(q2, temp, q); /* (q2*u) mod q */ 
mp_mult(temp, p, q2); /* p * [(q2*u) mod q] */ 
mp_add(expout, temp); /* expout = p2 + p * [...] */ 
#if 0 
} 
#endif 

mp_burn(q2); /* burn the evidence on the stack... */ 
mp_burn(temp); 
return 0; /* normal return */ 
} /* mp_modexp_crt */ 


/****************** end of MPI library ****************************/