/*
* Copyright (c) 1995 Colin Plumb. All rights reserved.
* For licensing and other legal details, see the file legal.c.
*
* sieve.c - Trial division for prime finding.
*
* Finding primes:
* - Sieve 1 to find the small primes for
* - Sieve 2 to find the candidate large primes, then
* - Pseudo-primality test.
*
* An important question is how much trial division by small primes
* should we do? The answer is a LOT. Even a heavily optimized
* Fermat test to the base 2 (the simplest pseudoprimality test)
* is much more expensive than a division.
*
* For an prime of n k-bit words, a Fermat test to the base 2 requires n*k
* modular squarings, each of which involves n*(n+1)/2 signle-word multiplies
* in the squaring and n*(n+1) multiplies in the modular reduction, plus
* some overhead to get into and out of Montgomery form. This is a total
* of 3/2 * k * n^2 * (n+1). Equivalently, if n*k = b bits, it's
* 3/2 * (b/k+1) * b^2 / k.
*
* A modulo operation requires n single-word divides. Let's assume that
* a divide is 4 times the cost of a multiply. That's 4*n multiplies.
* However, you only have to do the division once for your entire
* search. It can be amortized over 10-15 primes. So it's
* really more like n/3 multiplies. This is b/3k.
*
* Now, let's suppose you have a candidate prime t. Your options
* are to a) do trial division by a prime p, then do a Fermat test,
* or to do the Fermat test directly. Doing the trial division
* costs b/3k multiplies, but a certain fraction of the time (1/p), it
* saves you 3/2 b^3 / k^2 multiplies. Thus, it's worth it doing the
* division as long as b/3k < 3/2 * (b/k+1) * b^2 / k / p.
* I.e. p < 9/2 * (b/k + 1) * b = 9/2 * (b^2/k + b).
* E.g. for k=16 and b=256, p < 9/2 * 17 * 256 = 19584.
* Solving for k=16 and k=32 at a few interesting value of b:
*
* k=16, b=256: p < 19584 k=32, b=256: p < 10368
* k=16, b=384: p < 43200 k=32, b=384; p < 22464
* k=16, b=512: p < 76032 k=32, b=512: p < 39168
* k=16, b=640: p < 118080 k=32, b=640: p < 60480
*
* H'm... before using the highly-optimized Fermat test, I got much larger
* numbers (64K to 256K), and designed the sieve for that. Maybe it needs
* to be reduced. It *is* true that the desirable sieve size increases
* rapidly with increasing prime size, and it's the larger primes that are
* worrisome in any case. I'll leave it as is (64K) for now while I
* think about it.
*
* A bit of tweaking the division (we can compute a reciprocal and do
* multiplies instead, turning 4*n into 4 + 2*n) would increase all the
* numbers by a factor of 2 or so.
*
*
* Bit k in a sieve corresponds to the number a + k*b.
* For a given a and b, the sieve's job is to find the values of
* k for which a + k*b == 0 (mod p). Multiplying by b^-1 and
* isolating k, you get k == -a*b^-1 (mod p). So the values of
* k which should be worked on are k = (-a*b^-1 mod p) + i * p,
* for i = 0, 1, 2,...
*
* Note how this is still easy to use with very large b, if you need it.
* It just requires computing (b mod p) and then finding the multiplicative
* inverse of that.
*
*
* How large a space to search to ensure that one will hit a prime?
* The average density is known, but the primes behave oddly, and sometimes
* there are large gaps. It is conjectured by shanks that the first gap
* of size "delta" will occur at approximately exp(sqrt(delta)), so a delta
* of 65536 is conjectured to be to contain a prime up to e^256.
* Remembering the handy 2<->e conversion ratios:
* ln(2) = 0.693147 log2(e) = 1.442695
* This covers up to 369 bits. Damn, not enough! Still, it'll have to do.
*
* Cramer's conjecture (he proved it for "most" cases) is that in the limit,
* as p goes to infinity, the largest gap after a prime p tends to (ln(p))^2.
* So, for a 1024-bit p, the interval to the next prime is expected to be
* about 709.78^2, or 503791. We'd need to enlarge our space by a factor of
* 8 to be sure. It isn't worth the hassle.
*
* Note that a span of this size is expected to contain 92 primes even
* in the vicinity of 2^1024 (it's 369 at 256 bits and 492 at 192 bits).
* So the probability of failure is pretty low.
*/
#ifndef HAVE_CONFIG_H
#define HAVE_CONFIG_H 0
#endif
#if HAVE_CONFIG_H
#include "bnconfig.h"
#endif
/*
* Some compilers complain about #if FOO if FOO isn't defined,
* so do the ANSI-mandated thing explicitly...
*/
#ifndef NO_ASSERT_H
#define NO_ASSERT_H 0
#endif
#ifndef NO_LIMITS_H
#define NO_LIMITS_H 0
#endif
#ifndef NO_STRING_H
#define NO_STRING_H 0
#endif
#ifndef HAVE_STRINGS_H
#define HAVE_STRINGS_H 0
#endif
#if !NO_ASSERT_H
#include <assert.h>
#else
#define assert(x) (void)0
#endif
#if !NO_LIMITS_H
#include <limits.h> /* For UINT_MAX */
#endif /* If not avail, default value of 0 is safe */
#if !NO_STRING_H
#include <string.h> /* for memset() */
#elif HAVE_STRINGS_H
#include <strings.h>
#endif
#include "bn.h"
#include "sieve.h"
#ifdef MSDOS
#include "lbnmem.h"
#endif
#include "kludge.h"
/*
* Each array stores potential primes as 1 bits in little-endian bytes.
* Bit k in an array represents a + k*b, for some parameters a and b
* of the sieve. Currently, b is hardcoded to 2.
*
* Various factors of 16 arise because these are all *byte* sizes, and
* skipping even numbers, 16 numbers fit into a byte's worth of bitmap.
*/
/*
* The first number in the small prime sieve. This could be raised to
* 3 if you want to squeeze bytes out aggressively for a smaller SMALL
* table, and doing so would let one more prime into the end of the array,
* but there is no sense making it larger if you're generating small
* primes up to the limit if 2^16, since it doesn't save any memory and
* would require extra code to ignore 65537 in the last byte, which is
* over the 16-bit limit.
*/
#define SMALLSTART 1
/*
* Size of sieve used to find large primes, in bytes. For compatibility
* with 16-bit-int systems, the largest prime that can appear in it,
* SMALL * 16 + SMALLSTART - 2, must be < 65536. Since 65537 is a prime,
* this is the absolute maximum table size.
*/
#define SMALL (65536/16)
/*
* Compute the multiplicative inverse of x, modulo mod, using the extended
* Euclidean algorithm. The classical EEA returns two results, traditionally
* named s and t, but only one (t) is needed or computed here.
* It is unrolled twice to avoid some variable-swapping, and because negating
* t every other round makes all the number positive and less than the
* modulus, which makes fixed-length arithmetic easier.
*
* If gcd(x, mod) != 1, then this will return 0.
*/
static unsigned
sieveModInvert(unsigned x, unsigned mod)
{
unsigned y;
unsigned t0, t1;
unsigned q;
if (x <= 1)
return x; /* 0 and 1 are self-inverse */
/*
* The first round is simplified based on the
* initial conditions t0 = 1 and t1 = 0.
*/
t1 = mod / x;
y = mod % x;
if (y <= 1)
return y ? mod - t1 : 0;
t0 = 1;
do {
q = x / y;
x = x % y;
t0 += q * t1;
if (x <= 1)
return x ? t0 : 0;
q = y / x;
y = y % x;
t1 += q * t0;
} while (y > 1);
return y ? mod - t1 : 0;
}
/*
* Perform a single sieving operation on an array. Clear bits "start",
* "start+step", "start+2*step", etc. from the array, up to the size
* limit (in BYTES) "size". All of the arguments must fit into 16 bits
* for portability.
*
* This is the core of the sieving operation. In addition to being
* called from the sieving functions, it is useful to call directly if,
* say, you want to exclude primes congruent to 1 mod 3, or whatever.
* (Although in that case, it would be better to change the sieving to
* use a step size of 6 and start == 5 (mod 6).)
*
* Originally, this was inlined in the code below (with various checks
* turned off where they could be inferred from the environment), but it
* turns out that all the sieving is so fast that it makes a negligible
* speed difference and smaller, cleaner code was preferred.
*
* Rather than increment a bit index through the array and clear
* the corresponding bit, this code takes advantage of the fact that
* every eighth increment must use the same bit position in a byte.
* I.e. start + k*step == start + (k+8)*step (mod 8). Thus, a bitmask
* can be computed only eight times and used for all multiples. Thus, the
* outer loop is over (k mod 8) while the inner loop is over (k div 8).
*
* The only further trickiness is that this code is designed to accept
* start, step, and size up to 65535 on 16-bit machines. On such a
* machine, the computation "start+step" can overflow, so we need to
* insert an extra check for that situation.
*/
void
sieveSingle(unsigned char *array, unsigned size, unsigned start, unsigned step)
{
unsigned bit;
unsigned char mask;
unsigned i;
#if UINT_MAX < 0x1ffff
/* Unsigned is small; add checks for wrap */
for (bit = 0; bit < 8; bit++) {
i = start/8;
if (i >= size)
break;
mask = ~(1 << (start & 7));
do {
array[i] &= mask;
i += step;
} while (i >= step && i < size);
start += step;
if (start < step) /* Overflow test */
break;
}
#else
/* Unsigned has the range - no overflow possible */
for (bit = 0; bit < 8; bit++) {
i = start/8;
if (i >= size)
break;
mask = ~(1 << (start & 7));
do {
array[i] &= mask;
i += step;
} while (i < size);
start += step;
}
#endif
}
/*
* Returns the index of the next bit set in the given array. The search
* begins after the specified bit, so if you care about bit 0, you need
* to check it explicitly yourself. This returns 0 if no bits are found.
*
* Note that the size is in bytes, and that it takes and returns BIT
* positions. If the array represents odd numbers only, as usual, the
* returned values must be doubled to turn them into offsets from the
* initial number.
*/
unsigned
sieveSearch(unsigned char const *array, unsigned size, unsigned start)
{
unsigned i; /* Loop index */
unsigned char t; /* Temp */
if (!++start)
return 0;
i = start/8;
if (i >= size)
return 0; /* Done! */
/* Deal with odd-bit beginnings => search the first byte */
if (start & 7) {
t = array[i++] >> (start & 7);
if (t) {
if (!(t & 15)) {
t >>= 4;
start += 4;
}
if (!(t & 3)) {
t >>= 2;
start += 2;
}
if (!(t & 1))
start += 1;
return start;
} else if (i == size) {
return 0; /* Done */
}
}
/* Now the main search loop */
do {
if ((t = array[i]) != 0) {
start = 8*i;
if (!(t & 15)) {
t >>= 4;
start += 4;
}
if (!(t & 3)) {
t >>= 2;
start += 2;
}
if (!(t & 1))
start += 1;
return start;
}
} while (++i < size);
/* Failed */
return 0;
}
/*
* 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
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 */
}