bnsieve.c

著名的加密软件的应用于电子邮件中

/*
* Build a table of small primes for sieving larger primes with. This
* could be cached between calls to sieveBuild, but it's so fast that
* it's really not worth it. This code takes a few milliseconds to run.
*/
static void
 sieveSmall(unsigned char *array, unsigned size)
{
		unsigned i;		 /* Loop index */
		unsigned p;		 /* The current prime */

		/* Initialize to all 1s */
		memset(array, 0xFF, size);

#if SMALLSTART == 1
		/* Mark 1 as NOT prime */
		array[0] = 0xfe;
		i = 1;	/* Index of first prime */
#else
		i = 0;	/* Index of first prime */
#endif

		/*
		* Okay, now sieve via the primes up to 256, obtained from the
		* table itself. We know the maximum possible table size is
		* 65536, and sieveSingle() can cope with out-of-range inputs
		* safely, and the time required is trivial, so it isn't adaptive
		* based on the array size.
		*
		* Convert each bit position into a prime, compute a starting
		* sieve position (the square of the prime), and remove multiples
		* from the table, using sieveSingle(). I used to have that
		* code in line here, but the speed difference was so small it
		* wasn't worth it. If a compiler really wants to waste memory,
		* it can inline it.
		*/
		do {
				p = 2 * i + SMALLSTART;
				if (p > 256)
					break;
				/* Start at square of p */
				sieveSingle(array, size, (p*p-SMALLSTART)/2, p);

				/* And find the next prime */
				i = sieveSearch(array, 16, i);
		} while (i);
}


/*
* This is the primary sieving function. It fills in the array with
* a sieve (multiples of small primes removed) beginning at bn and
* proceeding in steps of "step".
*
* It generates a small array to get the primes to sieve by. It's
* generated on the fly - sieveSmall is fast enough to make that
* perfectly acceptable.
*
* The caller should take the array, walk it with sieveSearch, and
* apply a stronger primality test to the numbers that are returned.
*
* If the "dbl" flag non-zero (at least 1), this also sieves 2*bn+1, in
* steps of 2*step. If dbl is 2 or more, this also sieve 4*bn+3,
* in steps of 4*step, and so on for arbitrarily high values of "dbl".
* This is convenient for finding primes such that (p-1)/2 is also prime.
* This is particularly efficient because sieveSingle is controlled by the
* parameter s = -n/step (mod p). (In fact, we find t = -1/step (mod p)
* and multiply that by n (mod p).) If you have -n/step (mod p), then
* finding -(2*n+1)/(2*step) (mod p), which is -n/step - 1/(2*step) (mod p),
* reduces to finding -1/(2*step) (mod p), or t/2 (mod p), and adding that
* to s = -n/step (mod p). Dividing by 2 modulo an odd p is easy -
* if even, divide directly. Otherwise, add p (which produces an even
* sum), and divide by 2. Very simple. And this produces s' and t'
* for step' = 2*step. It can be repeated for step'' = 4*step and so on.
*
* Note that some of the math is complicated by the fact that 2*p might
* not fit into an unsigned, so rather than if (odd(x)) x = (x+p)/2,
* we do if (odd(x)) x = x/2 + p/2 + 1;
*
* TODO: Do the double-sieving by sieving the larger number, and then
* just subtract one from the remainder to get the other parameter.
* (bn-1)/2 is divisible by an odd p iff bn-1 is divisible, which is
* true iff bn == 1 mod p. This requires using a step size of 4.
*/
int
sieveBuild(unsigned char *array, unsigned size, struct BigNum const *bn,
		unsigned step, unsigned dbl)
{
		unsigned i, j;	/* Loop index */
		unsigned p;	/* Current small prime */
		unsigned s;	/* Where to start operations in the big sieve */
		unsigned t;	/* Step modulo p, the current prime */
#ifdef MSDOS	/* Use dynamic allocation rather than on the stack */
		unsigned char *small;
#else
		unsigned char small[SMALL];
#endif

		pgpAssert(array);

#ifdef MSDOS
		small = bniMemAlloc(SMALL);		/* Which allocator? Not secure. */
		if (!small)
			return -1;	/* Failed */
#endif

/*
* An odd step is a special case, since we must sieve by 2,
* which isn't in the small prime array and has a few other
* special properties. These are:
* - Since the numbers are stored in binary, we don't need to
* use bnModQ to find the remainder.
* - If step is odd, then t = step % 2 is 1, which allows
* the elimination of a lot of math. Inverting and negating
* t don't change it, and multiplying s by 1 is a no-op,
* so t isn't actually mentioned.
* - Since this is the first sieving, instead of calling
* sieveSingle, we can just use memset to fill the array
* with 0x55 or 0xAA. Since a 1 bit means possible prime
* (i.e. NOT divisible by 2), and the least significant bit
* is first, if bn % 2 == 0, we use 0xAA (bit 0 = bn is NOT
* prime), while if bn % 2 == 1, use 0x55.
* (If step is even, bn must be odd, so fill the array with 0xFF.)
* - Any doublings need not be considered, since 2*bn+1 is odd, and
* 2*step is even, so none of these numbers are divisible by 2.
*/
		if (step & 1) {
				s = bnLSWord(bn) & 1;
				memset(array, 0xAA >> s, size);
		} else {
				/* Initialize the array to all 1's */
				memset(array, 255, size);
				pgpAssert(bnLSWord(bn) & 1);
}

/*
* This could be cached between calls to sieveBuild, but
* it's really not worth it; sieveSmall is *very* fast.
* sieveSmall returns a sieve of odd primes.
*/
		sieveSmall(small, SMALL);

/*
* Okay, now sieve via the primes up to ssize*16+SMALLSTART-1,
* obtained from the small table.
*/
		i = (small[0] & 1) ? 0 : sieveSearch(small, SMALL, 0);
				do {
				p = 2 * i + SMALLSTART;

		/*
		* Modulo is usually very expensive, but step is usually
		* small, so this conditional is worth it.
		*/
		t = (step < p) ? step : step % p;
		if (!t) {
			/*
			* Instead of pgpAssert failing, returning all zero
			* bits is the "correct" thing to do, but I think
			* that the caller should take care of that
			* themselves before starting.
			*/
			pgpAssert(bnModQ(bn, p) != 0);
			continue;
	}
	/*
	* Get inverse of step mod p. 0 < t < p, and p is prime,
	* so it has an inverse and sieveModInvert can't return 0.
	*/
			t = sieveModInvert(t, p);
			pgpAssert(t);
			/* Negate t, so now t == -1/step (mod p) */
			t = p - t;

			/* Now get the bignum modulo the prime. */
			s = bnModQ(bn, p);

			/* Multiply by t, the negative inverse of step size */
#if UINT_MAX/0xffff < 0xffff
		s = (unsigned)(((unsigned long)s * t) % p);
#else
		s = (s * t) % p;
#endif

				/* s is now the starting bit position, so sieve */
				sieveSingle(array, size, s, p);

				/* Now do the double sieves as desired. */
				for (j = 0; j < dbl; j++) {
						/* Halve t modulo p */
#if UINT_MAX < 0x1ffff
						t = (t & 1) ? p/2 + t/2 + 1 : t/2;
						/* Add t to s, modulo p with overflow checks. */
						s += t;
						if (s >= p || s < t)
							s -= p;
#else
						if (t & 1)
							t += p;
						t /= 2;
						/* Add t to s, modulo p */
						s += t;
						if (s >= p)
							s -= p;
#endif
						sieveSingle(array, size, s, p);
				}

				/* And find the next prime */
		} while ((i = sieveSearch(small, SMALL, i)) != 0);

#ifdef MSDOS
		bniMemFree(small, SMALL);
#endif
		return 0;	 /* Success */
}

/*
* Similar to the above, but use "step" (which must be even) as a step
* size rather than a fixed value of 2. If "step" has any small divisors
* other than 2, this will blow up.
*
* Returns -1 on out of memory (MSDOS only, actually), and -2
* if step is found to be non-prime.
*/
int
sieveBuildBig(unsigned char *array, unsigned size, struct BigNum const *bn,
		struct BigNum const *step, unsigned dbl)
{
		unsigned i, j;	/* Loop index */
		unsigned p;	/* Current small prime */
		unsigned s;	/* Where to start operations in the big sieve */
		unsigned t;	/* step modulo p, the current prime */
#ifdef MSDOS	/* Use dynamic allocation rather than on the stack */
	unsigned char *small;
#else
	unsigned char small[SMALL];
#endif

	pgpAssert(array);

#ifdef MSDOS
		small = bniMemAlloc(SMALL);	/* Which allocator? Not secure. */
		if (!small)
			return -1;	/* Failed */
#endif
		/*
		* An odd step is a special case, since we must sieve by 2,
		* which isn't in the small prime array and has a few other
		* special properties. These are:
		* - Since the numbers are stored in binary, we don't need to
		* use bnModQ to find the remainder.
		* - If step is odd, then t = step % 2 is 1, which allows
		* the elimination of a lot of math. Inverting and negating
		* t don't change it, and multiplying s by 1 is a no-op,
		* so t isn't actually mentioned.
		* - Since this is the first sieving, instead of calling
		* sieveSingle, we can just use memset to fill the array
		* with 0x55 or 0xAA. Since a 1 bit means possible prime
		* (i.e. NOT divisible by 2), and the least significant bit
		* is first, if bn % 2 == 0, we use 0xAA (bit 0 = bn is NOT
		* prime), while if bn % 2 == 1, use 0x55.
		* (If step is even, bn must be odd, so fill the array with 0xFF.)
		* - Any doublings need not be considered, since 2*bn+1 is odd, and
		* 2*step is even, so none of these numbers are divisible by 2.
		*/
		if (bnLSWord(step) & 1) {
			s = bnLSWord(bn) & 1;
			memset(array, 0xAA >> s, size);
		} else {
				/* Initialize the array to all 1's */
				memset(array, 255, size);
				pgpAssert(bnLSWord(bn) & 1);
		}

/*
* This could be cached between calls to sieveBuild, but
* it's really not worth it; sieveSmall is *very* fast.
* sieveSmall returns a sieve of the odd primes.
*/
sieveSmall(small, SMALL);

/*
* Okay, now sieve via the primes up to ssize*16+SMALLSTART-1,
* obtained from the small table.
*/
		i = (small[0] & 1) ? 0 : sieveSearch(small, SMALL, 0);
		do {
				p = 2 * i + SMALLSTART;

				t = bnModQ(step, p);
				if (!t) {
					pgpAssert(bnModQ(bn, p) != 0);
					continue;
				}
				/* Get negative inverse of step */
				t = sieveModInvert(bnModQ(step, p), p);
				pgpAssert(t);
				t = p-t;

				/* Okay, we have a prime - get the remainder */
				s = bnModQ(bn, p);

				/* Now multiply s by the negative inverse of step (mod p) */
#if UINT_MAX < 0xffff * 0xffff
				s = (unsigned)(((unsigned long)s * t) % p);
#else
				s = (s * t) % p;
#endif
				/* We now have the starting bit pos */
				sieveSingle(array, size, s, p);

				/* Now do the double sieves as desired. */
				for (j = 0; j < dbl; j++) {
						/* Halve t modulo p */
#if UINT_MAX < 0x1ffff
						t = (t & 1) ? p/2 + t/2 + 1 : t/2;
						/* Add t to s, modulo p with overflow checks. */
						s += t;
						if (s >= p || s < t)
							s -= p;
#else
						if (t & 1)
							t += p;
						t /= 2;
						/* Add t to s, modulo p */
						s += t;
						if (s >= p)
							s -= p;
#endif
						sieveSingle(array, size, s, p);
				}

				/* And find the next prime */
		} while ((i = sieveSearch(small, SMALL, i)) != 0);

#ifdef MSDOS
	bniMemFree(small, SMALL);
#endif
	return 0;	 /* Success */
}