/*
 * sieve.c - Trial division for prime finding.
 *
 * Copyright (c) 1995  Colin Plumb.  All rights reserved.
 * For licensing and other legal details, see the file legal.c.
 *
 * 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
#ifndef NEED_MEMORY_H
#define NEED_MEMORY_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
#if NEED_MEMORY_H
#include <memory.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 */
}