genprime.cpp

RSA算法的VC源码

 * and p passes through the sieve, then p is definitely a prime.
 * If p is large and p passes through the sieve and may be a prime,
 * then p is further tested for primality with a slower test.
 */
{
    short i;
    static word16 lastprime = 0;	/* last prime in primetable */
	word16 sqrt_p;	/* to limit sieving past sqrt(p), for small p's */

    if (!lastprime) {		/* lastprime still undefined. So define it. */
	/* executes this code only once, then skips it next time */
	for (i = 0; primetable[i]; i++);	/* seek end of primetable */
	lastprime = primetable[i - 1];	/* get last prime in table */
    }
    if (significance(p) <= (2 / BYTES_PER_UNIT))	/* if p <= 16 bits */
	/* p may be in primetable.  Search it. */
	if (bottom16(p) <= lastprime)
	    for (i = 0; primetable[i]; i++) {
		/* scan until null-terminator */
		if (primetable[i] == bottom16(p))
			return TRUE;	/* yep, definitely a prime. */
		if (primetable[i] > bottom16(p))	/* we missed. */
			return FALSE;	/* definitely NOT a prime. */
	    }		/* We got past the whole primetable without a hit. */
    /* p is bigger than any prime in primetable, so let's sieve. */
    if (!(lsunit(p) & 1))	/* if least significant bit is 0... */
	return FALSE;		/* divisible by 2, not prime */

	if (mp_tstminus(p))		/* error if p<0 */
	return FALSE;		/* not prime if p<0 */

    /*
     * Optimization for small (32 bits or less) p's.  
     * If p is small, compute sqrt_p = sqrt(p), or else 
     * if p is >32 bits then just set sqrt_p to something 
     * at least as big as the largest primetable entry.
     */
    if (significance(p) <= (4 / BYTES_PER_UNIT)) {	/* if p <= 32 bits */
	unit sqrtp[MAX_UNIT_PRECISION];
	/* Just sieve up to sqrt(p) */
	if (mp_sqrt(sqrtp, p) == 0)	/* 0 means p is a perfect square */
	    return FALSE;	/* perfect square is not a prime */
	/* we know that sqrtp <= 16 bits because p <= 32 bits */
	sqrt_p = bottom16(sqrtp);
    } else {
	/* p > 32 bits, so obviate sqrt(p) test. */
	sqrt_p = lastprime;	/* ensures that we do ENTIRE sieve. */
	}

    /* p is assumed odd, so begin sieve at 3 */
    for (i = 1; primetable[i]; i++) {
	/* Compute p mod (primetable[i]).  If it divides evenly... */
	if (mp_shortmod(p, primetable[i]) == 0)
	    return FALSE;	/* then p is definitely NOT prime */
	if (primetable[i] > sqrt_p)	/* fully sieved p? */
	    return TRUE; /* yep, fully passed sieve, definitely a prime. */
    }
    /* It passed the sieve, so p is a suspected prime. */

    /*  Now try slow complex primality test on suspected prime. */
	return slowtest(p);		/* returns TRUE or FALSE */
}				/* primetest */

#endif

/*
 * Used in conjunction with fastsieve.  Builds a table of remainders 
 * relative to the random starting point p, so that fastsieve can
 * sequentially sieve for suspected primes quickly.  Call buildsieve
 * once, then call fastsieve for consecutive prime candidates.
 * Note that p must be odd, because the sieve begins at 3. 
 */
static void buildsieve(unitptr p, word16 remainders[])
{
    short i;
    for (i = 1; primetable[i]; i++) {
	remainders[i] = mp_shortmod(p, primetable[i]);
    }
}				/* buildsieve */

/*
   Fast prime sieving algorithm by Philip Zimmermann, March 1987.
   This is the fastest algorithm I know of for quickly sieving for
   large (hundreds of bits in length) "random" probable primes, because 
   it uses only single-precision (16-bit) arithmetic.  Because rigorous 
   prime testing algorithms are very slow, it is recommended that 
   potential prime candidates be quickly passed through this fast 
   sieving algorithm first to weed out numbers that are obviously not
   prime.

   This algorithm is optimized to search sequentially for a large prime 
   from a random starting point.  For generalized nonsequential prime 
   testing, the slower  conventional sieve should be used, as given 
   in primetest(p).

   This algorithm requires a fixed table (called primetable) of the 
   first hundred or so small prime numbers. 
   First we select a random odd starting point (p) for our prime 
   search.  Then we build a table of 16-bit remainders calculated
   from that initial p.  This table of 16-bit remainders is exactly 
   the same length as the table of small 16-bit primes.  Each
   remainders table entry contains the remainder of p divided by the 
   corresponding primetable entry.  Then we begin sequentially testing
   all odd integers, starting from the initial odd random p.  The 
   distance we have searched from the huge random starting point p is 
   a small 16-bit number, pdelta.  If pdelta plus a remainders table 
   entry is evenly divisible by the corresponding primetable entry, 
   then p+pdelta is factorable by that primetable entry, which means
   p+pdelta is not prime.
 */

/*      Fastsieve is used for searching sequentially from a random starting
   point for a suspected prime.  Requires that buildsieve be called 
   first, to build a table of remainders relative to the random starting 
   point p.  
   Returns true iff pdelta passes through the sieve and thus p+pdelta 
   may be a prime.  Note that p must be odd, because the sieve begins 
   at 3.
 */
static boolean fastsieve(word16 pdelta, word16 remainders[])//快速筛选
{
    short i;
	for (i = 1; primetable[i]; i++) {
	/*
	 * If pdelta plus a remainders table entry is evenly 
	 * divisible by the corresponding primetable entry,
	 * then p+pdelta is factorable by that primetable entry, 
	 * which means p+pdelta is not prime.
	 */
	if ((pdelta + remainders[i]) % primetable[i] == 0)
	    return FALSE;	/* then p+pdelta is not prime */
    }
    /* It passed the sieve.  It is now a suspected prime. */
    return TRUE;
}				/* fastsieve */


#define numberof(x) (sizeof(x)/sizeof(x[0]))	/* number of table entries */


static int nextprime(unitptr p)
/*
 * Find next higher prime starting at p, returning result in p.
 * Uses fast prime sieving algorithm to search sequentially.
 * Returns 0 for normal completion status, < 0 for failure status.
 */
{
	word16 pdelta, range;
	short oldprecision;
	short i, suspects;

	/* start search at candidate p */
	mp_inc(p);	/* remember, it's the NEXT prime from p, noninclusive. */
	if (significance(p) <= 1) {
	/*
	 * p might be smaller than the largest prime in primetable.
	 * We can't sieve for primes that are already in primetable.
	 * We will have to directly search the table.
	 */
	/* scan until null-terminator */
	for (i = 0; primetable[i]; i++) {
		if (primetable[i] >= lsunit(p)) {
		mp_init(p, primetable[i]);
		return 0;	/* return next higher prime from primetable */
		}
	}		/* We got past the whole primetable without a hit. */
	}	      /* p is bigger than any prime in primetable, so let's sieve. */
	if (mp_tstminus(p)) {	/* error if p<0 */
	mp_init(p, 2);		/* next prime >0 is 2 */
	return 0;		/* normal completion status */
	}
#ifndef BLUM
	lsunit(p) |= 1;		/* set candidate's lsb - make it odd */
#else
	lsunit(p) |= 3;		/* Make candidate ==3 mod 4 */
#endif

	/* 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(countbits(p) + 4 + SLOP_BITS));
	/* Rescale p to global_precision we just defined */
	rescale(p, oldprecision, global_precision);

	{
#ifdef _NOMALLOC /* No malloc and free functions available.  Use stack. */
	word16 remainders[numberof(primetable)];
#else			/* malloc available, we can conserve stack space. */
	word16 *remainders;
	/* Allocate some memory for the table of remainders: */
	remainders = (word16 *) malloc(sizeof(primetable));
#endif				/* malloc available */

	/* Build remainders table relative to initial p: */
	buildsieve(p, remainders);
	pdelta = 0;		/* offset from initial random p */
	/* Sieve preparation complete.  Now for some fast fast sieving... */
	/* slowtest will not be called unless fastsieve is true */

	/* range is how far to search before giving up. */
#ifndef BLUM
	range = 4 * units2bits(global_precision);
#else
	/* Twice as many because step size is twice as large, */
	range = 8 * units2bits(global_precision);
#endif
	suspects = 0;	/* number of suspected primes and slowtest trials */
	for (;;) {
		/* equivalent to:  if (primetest(p)) break; */
		if (fastsieve(pdelta, remainders)) { /* found suspected prime */
		suspects++;	/* tally for statistical purposes */
		if (slowtest(p))
			break;	/* found a prime */
		}
#ifndef BLUM
		pdelta += 2;	/* try next odd number */
#else
		pdelta += 4;
		mp_inc(p);
		mp_inc(p);
#endif
		mp_inc(p);
		mp_inc(p);
		if (pdelta > range)	/* searched too many candidates? */
		break;	/* something must be wrong--bail out of search */

	}			/* while (TRUE) */

	for (i = 0; primetable[i]; i++)	/* scan until null-terminator */
		remainders[i] = 0;	/* don't leave remainders exposed in RAM */
#ifndef _NOMALLOC
	free(remainders);	/* free allocated memory */
#endif				/* not _NOMALLOC */
	}

	set_precision(oldprecision);	/* restore precision */

	if (pdelta > range) {	/* searched too many candidates? */
	if (suspects < 1)	/* unreasonable to have found no suspects */
		return NOSUSPECTS;	/* fastsieve failed, probably */
	return NOPRIMEFOUND;	/* return error status */
	}
	return 0;			/* return normal completion status */

}				/* nextprime */


/* We will need a series of truly random bits for key generation.
   In most implementations, our random number supply is derived from
   random keyboard delays rather than a hardware random number
   chip.  So we will have to ensure we have a large enough pool of
   accumulated random numbers from the keyboard.  trueRandAccum()
   performs this operation.
 */

/* Fills 1 unit with random bytes, and returns unit. */
static unit randomunit(void)
{
	unit u = 0;
	byte i;
	i = BYTES_PER_UNIT;
	do
	u = (u << 8) +random(128); /* ZHAO YONG trueRandByte(); */
	while (--i != 0);
	return u;
}				/* randomunit */

/*
 * Make a random unit array p with nbits of precision.  Used mainly to 
 * generate large random numbers to search for primes.
 */
static void randombits(unitptr p, short nbits)
{
	//static int i=0;
	char r;
	mp_init(p, 0);
	make_lsbptr(p, global_precision);
	/* Add whole units of randomness */
	/*printf("Please Input some charater:");
	while (nbits >= UNITSIZE) {
			r=getchar();
			*post_higherunit(p) =r|randomunit();
			nbits -= UNITSIZE;
	}*/
	while(nbits >= UNITSIZE)
	{
		//r = InputChar[i++];
		r = rand()%0x80;//lzh
		*post_higherunit(p) =r|randomunit();
		nbits -= UNITSIZE;
	}
	/* Add most-significant partial unit (if any) */
	if (nbits)	*p = randomunit() & (power_of_2(nbits) - 1);
}				/* randombits */

/*
 * Makes a "random" prime p with nbits significant bits of precision.
 * Since these primes are used to compute a modulus of a guaranteed
 * length, the top 2 bits of the prime are set to 1, so that the
 * product of 2 primes (the modulus) is of a deterministic length.
 * Returns 0 for normal completion status, < 0 for failure status.
 */

int randomprime(unitptr p, short nbits)