genprime.cpp

RSA算法的VC源码

{
	DEBUGprintf2("\nGenerating a %d-bit random prime. ", nbits);
	/* Get an initial random candidate p to start search. */
	randombits(p, nbits - 2);	/* 2 less random bits for nonrandom top bits */
	/* To guarantee exactly nbits of significance, set the top 2 bits to 1 */
	/* ZHAO YONG ADD BEGIN */
	/*for( i=0;i<nbits/16;i++) p[i]=random(65535);
	p[nbits/16-1]=p[nbits/16-1]%16384;
	for(i=nbits/16;i<1024/16;i++) p[i]=0;
	*/
	// ZHAO YONG ADD END
	mp_setbit(p, nbits - 1);	/* highest bit is nonrandom */
	mp_setbit(p, nbits - 2);	/* next highest bit is also nonrandom */
	return nextprime(p);	/* search for next higher prime
				   from starting point p */
}				/* randomprime */


#ifdef STRONGPRIMES		/* generate "strong" primes for keys */

#define log_1stprime 6		/* log base 2 of firstprime */

/* 1st primetable entry used by tryprime */
#define firstprime (1<<log_1stprime)

/* This routine attempts to generate a prime p such that p-1 has prime p1
   as its largest factor.  Prime p will have no more than maxbits bits of
   significance.  Prime p1 must be less than maxbits-log_1stprime in length.  
   This routine is called only from goodprime.
 */
static boolean tryprime(unitptr p, unitptr p1, short maxbits)
{
    int i;
    unit i2[MAX_UNIT_PRECISION];
    /* Generate p such that p = (i*2*p1)+1, for i=1,2,3,5,7,11,13,17...
       and test p for primality for each small prime i.
       It's better to start i at firstprime rather than at 1,
	   because then p grows slower in significance.
       Start looking for small primes that are > firstprime...
     */
    if ((countbits(p1) + log_1stprime) >= maxbits) {
	DEBUGprintf1("\007[Error: overconstrained prime]");
	return FALSE;		/* failed to make a good prime */
    }
    for (i = 0; primetable[i]; i++) {
	if (primetable[i] < firstprime)
	    continue;
	/* note that mp_init doesn't extend sign bit for >32767 */
	mp_init(i2, primetable[i] << 1);
	mp_mult(p, p1, i2);
	mp_inc(p);
	if (countbits(p) > maxbits)
	    break;
	DEBUGprintf1(".");
	if (primetest(p))
	    return TRUE;
    }
	return FALSE;		/* failed to make a good prime */
}				/* tryprime */


/*
 * Make a "strong" prime p with at most maxbits and at least minbits of 
 * significant bits of precision.  This algorithm is called to generate
 * a high-quality prime p for key generation purposes.  It must have 
 * special characteristics for making a modulus n that is hard to factor.
 * Returns 0 for normal completion status, < 0 for failure status.
 */
int goodprime(unitptr p, short maxbits, short minbits)
{
    unit p1[MAX_UNIT_PRECISION];
	short oldprecision, midbits;
    int status;

    mp_init(p, 0);
    /* Adjust the global_precision downward to the optimum size for p... */
    oldprecision = global_precision;	/* save global_precision */
	/* We will need 2-3 extra bits of precision for the falsekeytest. */
    set_precision(bits2units(maxbits + 4 + SLOP_BITS));
    /* rescale p to global_precision we just defined */
    rescale(p, oldprecision, global_precision);

    minbits -= 2 * log_1stprime;	/* length of p" */
    midbits = (maxbits + minbits) / 2;	/* length of p' */
    DEBUGprintf3("\nGenerating a %d-%d bit refined prime. ",
		 minbits + 2 * log_1stprime, maxbits);
    do {
	do {
	    status = randomprime(p, minbits - 1);
		if (status < 0)
		return status;	/* failed to find a random prime */
		DEBUGprintf2("\n(p\042=%d bits)", countbits(p));
	} while (!tryprime(p1, p, midbits));
	DEBUGprintf2("(p'=%d bits)", countbits(p1));
    } while (!tryprime(p, p1, maxbits));
    DEBUGprintf2("\n\007(p=%d bits) ", countbits(p));
	set_precision(oldprecision);	/* restore precision */
    return 0;			/* normal completion status */
}				/* goodprime */

#endif				/* STRONGPRIMES */


#define iplus1  ( i==2 ? 0 : i+1 )	/* used by Euclid algorithms */
#define iminus1 ( i==0 ? 2 : i-1 )	/* used by Euclid algorithms */

/* Computes greatest common divisor via Euclid's algorithm. */
void mp_gcd(unitptr result, unitptr a, unitptr n)
{
    short i;
	unit gcopies[3][MAX_UNIT_PRECISION];
#define g(i) (  &(gcopies[i][0])  )
    mp_move(g(0), n);
    mp_move(g(1), a);

    i = 1;
	while (testne(g(i), 0)) {
	mp_mod(g(iplus1), g(iminus1), g(i));
	i = iplus1;
    }
    mp_move(result, g(iminus1));
#undef g
}				/* mp_gcd */

/*
 * Euclid's algorithm extended to compute multiplicative inverse.
 * Computes x such that a*x mod n = 1, where 0<a<n
 *
 * The variable u is unnecessary for the algorithm, but is
 * included in comments for mathematical clarity. 
 */
void mp_inv(unitptr x, unitptr a, unitptr n)
{
    short i;
	unit y[MAX_UNIT_PRECISION], temp[MAX_UNIT_PRECISION];
	unit gcopies[3][MAX_UNIT_PRECISION], vcopies[3][MAX_UNIT_PRECISION];
#define g(i) (  &(gcopies[i][0])  )
#define v(i) (  &(vcopies[i][0])  )
/*      unit ucopies[3][MAX_UNIT_PRECISION]; */
/* #define u(i) (  &(ucopies[i][0])  ) */
	mp_move(g(0), n);
	mp_move(g(1), a);
/*      mp_init(u(0),1); mp_init(u(1),0); */
	mp_init(v(0), 0);
	mp_init(v(1), 1);
	i = 1;
	while (testne(g(i), 0)) {
		/* we know that at this point,  g(i) = u(i)*n + v(i)*a  */
		mp_udiv(g(iplus1), y, g(iminus1), g(i));
		mp_mult(temp, y, v(i));
		mp_move(v(iplus1), v(iminus1));
		mp_sub(v(iplus1), temp);
		/*      mp_mult(temp,y,u(i)); mp_move(u(iplus1),u(iminus1));
		mp_sub(u(iplus1),temp); */
		i = iplus1;
	}
	mp_move(x, v(iminus1));
	if (mp_tstminus(x))	mp_add(x, n);
#undef g
#undef v
}				/* mp_inv */

#ifdef STRONGPRIMES

/*      mp_sqrt - returns square root of a number.
   returns -1 for error, 0 for perfect square, 1 for not perfect square.
   Not used by any RSA-related functions.       Used by factoring algorithms.
   This version needs optimization.
   by Charles W. Merritt  July 15, 1989, refined by PRZ.

   These are notes on computing the square root the manual old-fashioned 
   way.  This is the basis of the fast sqrt algorithm mp_sqrt below:

   1)   Separate the number into groups (periods) of two digits each,
   beginning with units or at the decimal point.
   2)   Find the greatest perfect square in the left hand period & write 
   its  square root as the first figure of the required root.  Subtract
   the square of this number from the left hand period.  Annex to the
   remainder the next group so as to form a dividend.
   3)   Double the root already found and write it as a partial divisor at 
   the left of the new dividend.  Annex one zero digit, making a trial
   divisor, and divide the new dividend by the trial divisor.
   4)   Write the quotient in the root as the trial term and also substitute 
   this quotient for the annexed zero digit in the partial divisor, 
   making the latter complete.
   5)   Multiply the complete divisor by the figure just obtained and, 
   if possible, subtract the product from the last remainder.
   If this product is too large, the trial term of the quotient
   must be replaced by the next smaller number and the operations
   preformed as before.
   (IN BINARY, OUR TRIAL TERM IS ALWAYS 1 AND WE USE IT OR DO NOTHING.)
   6)   Proceed in this manner until all periods are used.
   If there is still a remainder, it's not a perfect square.
 */

/* Quotient is returned as the square root of dividend. */
static int mp_sqrt(unitptr quotient, unitptr dividend)
{
	register short next2bits;	/* "period", or group of 2 bits of dividend */
    register unit dvdbitmask, qbitmask;
	unit remainder[MAX_UNIT_PRECISION];
    unit rjq[MAX_UNIT_PRECISION], divisor[MAX_UNIT_PRECISION];
    unsigned int qbits, qprec, dvdbits, dprec, oldprecision;
    int notperfect;

    mp_init(quotient, 0);
	if (mp_tstminus(dividend)) {	/* if dividend<0, return error */
	mp_dec(quotient);	/* quotient = -1 */
	return -1;
    }
    /* normalize and compute number of bits in dividend first */
    init_bitsniffer(dividend, dvdbitmask, dprec, dvdbits);
    /* init_bitsniffer returns a 0 if dvdbits is 0 */
    if (dvdbits == 1) {
	mp_init(quotient, 1);	/* square root of 1 is 1 */
	return 0;
    }
    /* rescale quotient to half the precision of dividend */
	qbits = (dvdbits + 1) >> 1;
    qprec = bits2units(qbits);
	rescale(quotient, global_precision, qprec);
    make_msbptr(quotient, qprec);
    qbitmask = power_of_2((qbits - 1) & (UNITSIZE - 1));

    /*
     * Set smallest optimum precision for this square root.
	 * The low-level primitives are affected by the call to set_precision.
     * Even though the dividend precision is bigger than the precision
     * we will use, no low-level primitives will be used on the dividend.
     * They will be used on the quotient, remainder, and rjq, which are
     * smaller precision.
     */
    oldprecision = global_precision;	/* save global_precision */
    set_precision(bits2units(qbits + 3));	/* 3 bits of precision slop */

    /* special case: sqrt of 1st 2 (binary) digits of dividend
       is 1st (binary) digit of quotient.  This is always 1. */
    stuff_bit(quotient, qbitmask);
	bump_bitsniffer(quotient, qbitmask);
    mp_init(rjq, 1);		/* rjq is Right Justified Quotient */

    if (!(dvdbits & 1)) {
	/* even number of bits in dividend */
	next2bits = 2;
	bump_bitsniffer(dividend, dvdbitmask);
	dvdbits--;
	if (sniff_bit(dividend, dvdbitmask))
	    next2bits++;
	bump_bitsniffer(dividend, dvdbitmask);
	dvdbits--;
    } else {
	/* odd number of bits in dividend */
	next2bits = 1;
	bump_bitsniffer(dividend, dvdbitmask);
	dvdbits--;
    }

    mp_init(remainder, next2bits - 1);

    /* dvdbits is guaranteed to be even at this point */

    while (dvdbits) {
	next2bits = 0;
	if (sniff_bit(dividend, dvdbitmask))
	    next2bits = 2;
	bump_bitsniffer(dividend, dvdbitmask);
	dvdbits--;
	if (sniff_bit(dividend, dvdbitmask))
	    next2bits++;
	bump_bitsniffer(dividend, dvdbitmask);
	dvdbits--;
	mp_rotate_left(remainder, (boolean) ((next2bits & 2) != 0));
	mp_rotate_left(remainder, (boolean) ((next2bits & 1) != 0));

	/*
	 * "divisor" is trial divisor, complete divisor is 4*rjq 
	 * or 4*rjq+1.
	 * Subtract divisor times its last digit from remainder.
	 * If divisor ends in 1, remainder -= divisor*1,
	 * or if divisor ends in 0, remainder -= divisor*0 (do nothing).
	 * Last digit of divisor inflates divisor as large as possible
	 * yet still subtractable from remainder.
	 */
	mp_move(divisor, rjq);	/* divisor = 4*rjq+1 */
	mp_rotate_left(divisor, 0);
	mp_rotate_left(divisor, 1);
	if (mp_compare(remainder, divisor) >= 0) {
	    mp_sub(remainder, divisor);
	    stuff_bit(quotient, qbitmask);
	    mp_rotate_left(rjq, 1);
	} else {
	    mp_rotate_left(rjq, 0);
	}
	bump_bitsniffer(quotient, qbitmask);
    }
    notperfect = testne(remainder, 0);	/* not a perfect square? */
    set_precision(oldprecision);	/* restore original precision */
    return notperfect;		/* normal return */
}				/* mp_sqrt */
#endif