sieve.c

包含哈希,对称以及非对称的经典算法 包含经典事例

	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

	assert(array);

#ifdef MSDOS
	small = lbnMemAlloc(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);
		assert(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 assert 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.
			 */
			assert(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);
		assert(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
	lbnMemFree(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

	assert(array);

#ifdef MSDOS
	small = lbnMemAlloc(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);
		assert(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) {
			assert(bnModQ(bn, p) != 0);
			continue;
		}
		/* Get negative inverse of step */
		t = sieveModInvert(bnModQ(step, p), p);
		assert(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
	lbnMemFree(small, SMALL);
#endif
	return 0;	/* Success */
}