mpilib.cpp

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

    mutemp = (MULTUNIT) min(q, (unsigned long)(1L << MULTUNITSIZE) - 1);
    q = q0 = q1 = q2 = 0;
    lsb_factor = 0;
    (void) lsb_factor;
    return mutemp;
}
#endif				/* mp_quo_digit */

/*
 * stage_smith_modulus() - Prepare for a Smith modmult.
 * 
 * Calculate the reciprocal of modulus with a precision of two MULTUNITs.
 * Assumes that global_precision has already been adjusted to the
 * size of the modulus, plus SLOP_BITS.
 * 
 * Note: This routine was designed to work with large values and
 * doesn't have the necessary testing or handling to work with a
 * modulus having less than three significant units.  For such cases,
 * the separate multiply and modulus routines can be used.
 * 
 * stage_smith_modulus() is aliased to stage_modulus().
 */
int stage_smith_modulus(const unit * n_modulus)
{
    int original_precision;
    int sigmod;			/* significant units in modulus */
    unitptr mp;			/* modulus most significant pointer */
    MULTUNIT *mpm;		/* reciprocal pointer */
    int prec;			/* precision of reciprocal calc in units */

    mp_move(modulus, n_modulus);
    modmpl = (MULTUNIT *) modulus;
    modmpl = lsbptr(modmpl, global_precision * MULTUNITs_perunit);
    nbits = countbits(modulus);
    modlenMULTUNITS = (nbits + MULTUNITSIZE - 1) / MULTUNITSIZE;

    original_precision = global_precision;

    /* The following code copies the three most significant units of
     * modulus to mod_divisor.
     */
    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 & 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_i