home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / crypl200.zip / BNLIB / SIEVE.C < prev    next >
C/C++ Source or Header  |  1996-05-16  |  21KB  |  689 lines
  1. /*
  2.  * sieve.c - Trial division for prime finding.
  3.  *
  4.  * Copyright (c) 1995  Colin Plumb.  All rights reserved.
  5.  * For licensing and other legal details, see the file legal.c.
  6.  *
  7.  * Finding primes:
  8.  * - Sieve 1 to find the small primes for
  9.  * - Sieve 2 to find the candidate large primes, then
  10.  * - Pseudo-primality test.
  11.  *
  12.  * An important question is how much trial division by small primes
  13.  * should we do?  The answer is a LOT.  Even a heavily optimized
  14.  * Fermat test to the base 2 (the simplest pseudoprimality test)
  15.  * is much more expensive than a division.
  16.  *
  17.  * For an prime of n k-bit words, a Fermat test to the base 2 requires n*k
  18.  * modular squarings, each of which involves n*(n+1)/2 signle-word multiplies
  19.  * in the squaring and n*(n+1) multiplies in the modular reduction, plus
  20.  * some overhead to get into and out of Montgomery form.  This is a total
  21.  * of 3/2 * k * n^2 * (n+1).  Equivalently, if n*k = b bits, it's
  22.  * 3/2 * (b/k+1) * b^2 / k.
  23.  *
  24.  * A modulo operation requires n single-word divides.  Let's assume that
  25.  * a divide is 4 times the cost of a multiply.  That's 4*n multiplies.
  26.  * However, you only have to do the division once for your entire
  27.  * search.  It can be amortized over 10-15 primes.  So it's
  28.  * really more like n/3 multiplies.  This is b/3k.
  29.  *
  30.  * Now, let's suppose you have a candidate prime t.  Your options
  31.  * are to a) do trial division by a prime p, then do a Fermat test,
  32.  * or to do the Fermat test directly.  Doing the trial division
  33.  * costs b/3k multiplies, but a certain fraction of the time (1/p), it
  34.  * saves you 3/2 b^3 / k^2 multiplies.  Thus, it's worth it doing the
  35.  * division as long as b/3k < 3/2 * (b/k+1) * b^2 / k / p.
  36.  * I.e. p < 9/2 * (b/k + 1) * b = 9/2 * (b^2/k + b).
  37.  * E.g. for k=16 and b=256, p < 9/2 * 17 * 256 = 19584.
  38.  * Solving for k=16 and k=32 at a few interesting value of b:
  39.  *
  40.  * k=16, b=256: p <  19584    k=32, b=256: p <  10368
  41.  * k=16, b=384: p <  43200    k=32, b=384; p <  22464
  42.  * k=16, b=512: p <  76032    k=32, b=512: p <  39168
  43.  * k=16, b=640: p < 118080    k=32, b=640: p <  60480
  44.  *
  45.  * H'm... before using the highly-optimized Fermat test, I got much larger
  46.  * numbers (64K to 256K), and designed the sieve for that.  Maybe it needs
  47.  * to be reduced.  It *is* true that the desirable sieve size increases
  48.  * rapidly with increasing prime size, and it's the larger primes that are
  49.  * worrisome in any case.  I'll leave it as is (64K) for now while I
  50.  * think about it.
  51.  *
  52.  * A bit of tweaking the division (we can compute a reciprocal and do
  53.  * multiplies instead, turning 4*n into 4 + 2*n) would increase all the
  54.  * numbers by a factor of 2 or so.
  55.  *
  56.  *
  57.  * Bit k in a sieve corresponds to the number a + k*b.
  58.  * For a given a and b, the sieve's job is to find the values of
  59.  * k for which a + k*b == 0 (mod p).  Multiplying by b^-1 and
  60.  * isolating k, you get k == -a*b^-1 (mod p).  So the values of
  61.  * k which should be worked on are k = (-a*b^-1 mod p) + i * p,
  62.  * for i = 0, 1, 2,...
  63.  *
  64.  * Note how this is still easy to use with very large b, if you need it.
  65.  * It just requires computing (b mod p) and then finding the multiplicative
  66.  * inverse of that.
  67.  *
  68.  *
  69.  * How large a space to search to ensure that one will hit a prime?
  70.  * The average density is known, but the primes behave oddly, and sometimes
  71.  * there are large gaps.  It is conjectured by shanks that the first gap
  72.  * of size "delta" will occur at approximately exp(sqrt(delta)), so a delta
  73.  * of 65536 is conjectured to be to contain a prime up to e^256.
  74.  * Remembering the handy 2<->e conversion ratios:
  75.  * ln(2) = 0.693147   log2(e) = 1.442695
  76.  * This covers up to 369 bits.  Damn, not enough!  Still, it'll have to do.
  77.  *
  78.  * Cramer's conjecture (he proved it for "most" cases) is that in the limit,
  79.  * as p goes to infinity, the largest gap after a prime p tends to (ln(p))^2.
  80.  * So, for a 1024-bit p, the interval to the next prime is expected to be
  81.  * about 709.78^2, or 503791.  We'd need to enlarge our space by a factor of
  82.  * 8 to be sure.  It isn't worth the hassle.
  83.  *
  84.  * Note that a span of this size is expected to contain 92 primes even
  85.  * in the vicinity of 2^1024 (it's 369 at 256 bits and 492 at 192 bits).
  86.  * So the probability of failure is pretty low.
  87.  */
  88. #ifndef HAVE_CONFIG_H
  89. #define HAVE_CONFIG_H 0
  90. #endif
  91. #if HAVE_CONFIG_H
  92. #include "config.h"
  93. #endif
  94.  
  95. /*
  96.  * Some compilers complain about #if FOO if FOO isn't defined,
  97.  * so do the ANSI-mandated thing explicitly...
  98.  */
  99. #ifndef NO_ASSERT_H
  100. #define NO_ASSERT_H 0
  101. #endif
  102. #ifndef NO_LIMITS_H
  103. #define NO_LIMITS_H 0
  104. #endif
  105. #ifndef NO_STRING_H
  106. #define NO_STRING_H 0
  107. #endif
  108. #ifndef HAVE_STRINGS_H
  109. #define HAVE_STRINGS_H 0
  110. #endif
  111. #ifndef NEED_MEMORY_H
  112. #define NEED_MEMORY_H 0
  113. #endif
  114.  
  115. #if !NO_ASSERT_H
  116. #include <assert.h>
  117. #else
  118. #define assert(x) (void)0
  119. #endif
  120.  
  121. #if !NO_LIMITS_H
  122. #include <limits.h>    /* For UINT_MAX */
  123. #endif            /* If not avail, default value of 0 is safe */
  124.  
  125. #if !NO_STRING_H
  126. #include <string.h>    /* for memset() */
  127. #elif HAVE_STRINGS_H
  128. #include <strings.h>
  129. #endif
  130. #if NEED_MEMORY_H
  131. #include <memory.h>
  132. #endif
  133.  
  134. #include "bn.h"
  135. #include "sieve.h"
  136. #ifdef MSDOS
  137. #include "lbnmem.h"
  138. #endif
  139.  
  140. #include "kludge.h"
  141.  
  142. /*
  143.  * Each array stores potential primes as 1 bits in little-endian bytes.
  144.  * Bit k in an array represents a + k*b, for some parameters a and b
  145.  * of the sieve.  Currently, b is hardcoded to 2.
  146.  *
  147.  * Various factors of 16 arise because these are all *byte* sizes, and
  148.  * skipping even numbers, 16 numbers fit into a byte's worth of bitmap.
  149.  */
  150.  
  151. /*
  152.  * The first number in the small prime sieve.  This could be raised to 3
  153.  * if you want to squeeze bytes bytes out aggressively for a smaller SMALL
  154.  * table, and doing so would let one more prime into the end of the array,
  155.  * but there is no sense making it larger if you're generating small primes
  156.  * up to the limit if 2^16, since it doesn't save any memory and would
  157.  * require extra code to ignore 65537 in the last byte, which is over the
  158.  * 16-bit limit.
  159.  */
  160. #define SMALLSTART 1
  161.  
  162. /*
  163.  * Size of sieve used to find large primes, in bytes.  For compatibility
  164.  * with 16-bit-int systems, the largest prime that can appear in it,
  165.  * SMALL * 16 + SMALLSTART - 2, must be < 65536.  Since 65537 is a prime,
  166.  * this is the absolute maximum table size.
  167.  */
  168. #define SMALL (65536/16)
  169.  
  170. /*
  171.  * Compute the multiplicative inverse of x, modulo mod, using the extended
  172.  * Euclidean algorithm.  The classical EEA returns two results, traditionally
  173.  * named s and t, but only one (t) is needed or computed here.
  174.  * It is unrolled twice to avoid some variable-swapping, and because negating
  175.  * t every other round makes all the number positive and less than the
  176.  * modulus, which makes fixed-length arithmetic easier.
  177.  *
  178.  * If gcd(x, mod) != 1, then this will return 0.
  179.  */
  180. static unsigned
  181. sieveModInvert(unsigned x, unsigned mod)
  182. {
  183.     unsigned y;
  184.     unsigned t0, t1;
  185.     unsigned q;
  186.  
  187.     if (x <= 1)
  188.         return x;    /* 0 and 1 are self-inverse */
  189.     /*
  190.      * The first round is simplified based on the
  191.      * initial conditions t0 = 1 and t1 = 0.
  192.      */
  193.     t1 = mod / x;
  194.     y = mod % x;
  195.     if (y <= 1)
  196.         return y ? mod - t1 : 0;
  197.     t0 = 1;
  198.  
  199.     do {
  200.         q = x / y;
  201.         x = x % y;
  202.         t0 += q * t1;
  203.         if (x <= 1)
  204.             return x ? t0 : 0;
  205.         q = y / x;
  206.         y = y % x;
  207.         t1 += q * t0;
  208.     } while (y > 1);
  209.     return y ? mod - t1 : 0;
  210. }
  211.  
  212.  
  213. /*
  214.  * Perform a single sieving operation on an array.  Clear bits
  215.  * "start", "start+step", "start+2*step", etc. from the array,
  216.  * up to the size limit (in BYTES) "size".  All of the arguments
  217.  * must fit into 16 bits for portability.
  218.  *
  219.  * This is the core of the sieving operation.  In addition to being
  220.  * called from the sieving functions, it is useful to call directly
  221.  * if, say, you want to exclude primes congruent to 1 mod 3, or
  222.  * whatever.  (Although in that case, it would be better to change
  223.  * the sieving to use a step size of 6 and start == 5 (mod 6).)
  224.  *
  225.  * Originally, this was inlined in the code below (with various checks
  226.  * turned off where they could be inferred from the environment), but
  227.  * it turns out that all the sieving is so fast that it makes a
  228.  * negligible speed difference and smaller, cleaner code was preferred.
  229.  *
  230.  * Rather than increment a bit index through the array and clear
  231.  * the corresponding bit, this code takes advantage of the fact that
  232.  * every eighth increment must use the same bit position in a byte.
  233.  * I.e. start + k*step == start + (k+8)*step (mod 8).  Thus, a
  234.  * bitmask can be computed only eight times and used for all multiples.
  235.  * Thus, the outer loop is over (k mod 8) while the inner loop is
  236.  * over (k div 8).
  237.  *
  238.  * The only further trickiness is that this code is designed to accept
  239.  * start, step, and size up to 65535 on 16-bit machines.  On such
  240.  * a machine, the computation "start+step" can overflow, so we need
  241.  * to insert an extra check for that situation.
  242.  */
  243. void
  244. sieveSingle(unsigned char *array, unsigned size, unsigned start, unsigned step)
  245. {
  246.     unsigned bit;
  247.     unsigned char mask;
  248.     unsigned i;
  249.  
  250. #if UINT_MAX < 0x1ffff
  251.     /* Unsigned is small; add checks for wrap */
  252.     for (bit = 0; bit < 8; bit++) {
  253.         i = start/8;
  254.         if (i >= size)
  255.             break;
  256.         mask = ~(1 << (start & 7));
  257.         do {
  258.             array[i] &= mask;
  259.             i += step;
  260.         } while (i >= step && i < size);
  261.         start += step;
  262.         if (start < step)    /* Overflow test */
  263.             break;
  264.     }
  265. #else
  266.     /* Unsigned has the range - no overflow possible */
  267.     for (bit = 0; bit < 8; bit++) {
  268.         i = start/8;
  269.         if (i >= size)
  270.             break;
  271.         mask = ~(1 << (start & 7));
  272.         do {
  273.             array[i] &= mask;
  274.             i += step;
  275.         } while (i < size);
  276.         start += step;
  277.     }
  278. #endif
  279. }
  280.  
  281. /*
  282.  * Returns the index of the next bit set in the given array.
  283.  * The search begins after the specified bit, so if you care
  284.  * about bit 0, you need to check it explicitly yourself.
  285.  * This returns 0 if no bits are found.
  286.  *
  287.  * Note that the size is in bytes, and that it takes and returns
  288.  * BIT positions.  If the array represents odd numbers only,
  289.  * as usual, the returns values must be doubled to turn them into
  290.  * offsets from the initial number.
  291.  */
  292. unsigned
  293. sieveSearch(unsigned char const *array, unsigned size, unsigned start)
  294. {
  295.     unsigned i;    /* Loop index */
  296.     unsigned char t;    /* Temp */
  297.  
  298.     if (!++start)
  299.         return 0;
  300.     i = start/8;
  301.     if (i >= size)
  302.         return 0;    /* Done! */
  303.  
  304.     /* Deal with odd-bit beginnings => search the first byte */
  305.     if (start & 7) {
  306.         t = array[i++] >> (start & 7);
  307.         if (t) {
  308.             if (!(t & 15)) {
  309.                 t >>= 4;
  310.                 start += 4;
  311.             }
  312.             if (!(t & 3)) {
  313.                 t >>= 2;
  314.                 start += 2;
  315.             }
  316.             if (!(t & 1))
  317.                 start += 1;
  318.             return start;
  319.         } else if (i == size) {
  320.             return 0;    /* Done */
  321.         }
  322.     }
  323.  
  324.     /* Now the main search loop */
  325.  
  326.     do {
  327.         if ((t = array[i]) != 0) {
  328.             start = 8*i;
  329.             if (!(t & 15)) {
  330.                 t >>= 4;
  331.                 start += 4;
  332.             }
  333.             if (!(t & 3)) {
  334.                 t >>= 2;
  335.                 start += 2;
  336.             }
  337.             if (!(t & 1))
  338.                 start += 1;
  339.             return start;
  340.         }
  341.     } while (++i < size);
  342.  
  343.     /* Failed */
  344.     return 0;
  345. }
  346.  
  347. /*
  348.  * Build a table of small primes for sieving larger primes with.
  349.  * This could be cached between calls to sieveBuild, but it's so
  350.  * fast that it's really not worth it.  This code takes a few
  351.  * milliseconds to run.
  352.  */
  353. static void
  354. sieveSmall(unsigned char *array, unsigned size)
  355. {
  356.     unsigned i;        /* Loop index */
  357.     unsigned p;        /* The current prime */
  358.  
  359.     /* Initialize to all 1s */
  360.     memset(array, 0xFF, size);
  361.  
  362. #if SMALLSTART == 1
  363.     /* Mark 1 as NOT prime */
  364.     array[0] = 0xfe;
  365.     i = 1;    /* Index of first prime */
  366. #else
  367.     i = 0;    /* Index of first prime */
  368. #endif
  369.  
  370.     /*
  371.      * Okay, now sieve via the primes up to 256, obtained from the
  372.      * table itself.  We know the maximum possible table size
  373.      * is 65536, and sieveSingle() can cope with out-of-range inputs
  374.      * safely, and the time required is trivial, so it isn't adaptive
  375.      * based on the array size.
  376.      *
  377.      * Convert each bit position into a prime, compute a starting
  378.      * sieve position (the square of the prime), and remove multiples
  379.      * from the table, using sieveSingle().  I used to have that code
  380.      * in line here, but the speed difference was so small it wasn't
  381.      * worth it.  If a compiler really wants to waste memory, it
  382.      * can inline it.
  383.      */
  384.     do {
  385.         p = 2 * i + SMALLSTART;
  386.         if (p > 256)
  387.             break;
  388.         /* Start at square of p */
  389.         sieveSingle(array, size, (p*p-SMALLSTART)/2, p);
  390.  
  391.         /* And find the next prime */
  392.         i = sieveSearch(array, 16, i);
  393.     } while (i);
  394. }
  395.  
  396.  
  397. /*
  398.  * This is the primary sieving function.  It fills in the array with
  399.  * a sieve (multiples of small primes removed) beginning at bn and
  400.  * proceeding in steps of "step".
  401.  *
  402.  * It generates a small array to get the primes to sieve by.  It's
  403.  * generated on the fly - sieveSmall is fast enough to make that
  404.  * perfectly acceptable.
  405.  *
  406.  * The caller should take the array, walk it with sieveSearch, and
  407.  * apply a stronger primality test to the numbers that are returned.
  408.  *
  409.  * If the "dbl" flag non-zero (at least 1), this also sieves 2*bn+1, in
  410.  * steps of 2*step.  If dbl is 2 or more, this also sieve 4*bn+3,
  411.  * in steps of 4*step, and so on for arbitrarily high values of "dbl".
  412.  * This is convenient for finding primes such that (p-1)/2 is also prime.
  413.  * This is particularly efficient because sieveSingle is controlled by the
  414.  * parameter s = -n/step (mod p).  (In fact, we find t = -1/step (mod p)
  415.  * and multiply that by n (mod p).)  If you have -n/step (mod p), then
  416.  * finding -(2*n+1)/(2*step) (mod p), which is -n/step - 1/(2*step) (mod p),
  417.  * reduces to finding -1/(2*step) (mod p), or t/2 (mod p), and adding that
  418.  * to s = -n/step (mod p).  Dividing by 2 modulo an odd p is easy -
  419.  * if even, divide directly.  Otherwise, add p (which produces an even
  420.  * sum), and divide by 2.  Very simple.  And this produces s' and t'
  421.  * for step' = 2*step.  It can be repeated for step'' = 4*step and so on.
  422.  *
  423.  * Note that some of the math is complicated by the fact that 2*p might
  424.  * not fit into an unsigned, so rather than if (odd(x)) x = (x+p)/2,
  425.  * we do if (odd(x)) x = x/2 + p/2 + 1;
  426.  *
  427.  * TODO: Do the double-sieving by sieving the larger number, and then
  428.  * just subtract one from the remainder to get the other parameter.
  429.  * (bn-1)/2 is divisible by an odd p iff bn-1 is divisible, which is true
  430.  * iff bn == 1 mod p.  This requires using a step size of 4.
  431.  */
  432. int
  433. sieveBuild(unsigned char *array, unsigned size, struct BigNum const *bn,
  434.     unsigned step, unsigned dbl)
  435. {
  436.     unsigned i, j;    /* Loop index */
  437.     unsigned p;    /* Current small prime */
  438.     unsigned s;    /* Where to start operations in the big sieve */
  439.     unsigned t;    /* Step modulo p, the current prime */
  440. #ifdef MSDOS    /* Use dynamic allocation rather than on the stack */
  441.     unsigned char *small;
  442. #else
  443.     unsigned char small[SMALL];
  444. #endif
  445.  
  446.     assert(array);
  447.  
  448. #ifdef MSDOS
  449.     small = lbnMemAlloc(SMALL);    /* Which allocator?  Not secure. */
  450.     if (!small)
  451.         return -1;    /* Failed */
  452. #endif
  453.  
  454.     /*
  455.      * An odd step is a special case, since we must sieve by 2,
  456.      * which isn't in the small prime array and has a few other
  457.      * special properties.  These are:
  458.      * - Since the numbers are stored in binary, we don't need to
  459.      *   use bnModQ to find the remainder.
  460.      * - If step is odd, then t = step % 2 is 1, which allows
  461.      *   the elimination of a lot of math.  Inverting and negating
  462.      *   t don't change it, and multiplying s by 1 is a no-op,
  463.      *   so t isn't actually mentioned.
  464.      * - Since this is the first sieving, instead of calling
  465.      *   sieveSingle, we can just use memset to fill the array
  466.      *   with 0x55 or 0xAA.  Since a 1 bit means possible prime
  467.      *   (i.e. NOT divisible by 2), and the least significant bit
  468.      *   is first, if bn % 2 == 0, we use 0xAA (bit 0 = bn is NOT
  469.      *   prime), while if bn % 2 == 1, use 0x55.
  470.      *   (If step is even, bn must be odd, so fill the array with 0xFF.)
  471.      * - Any doublings need not be considered, since 2*bn+1 is odd, and
  472.      *   2*step is even, so none of these numbers are divisible by 2.
  473.      */
  474.     if (step & 1) {
  475.         s = bnLSWord(bn) & 1;
  476.         memset(array, 0xAA >> s, size);
  477.     } else {
  478.         /* Initialize the array to all 1's */
  479.         memset(array, 255, size);
  480.         assert(bnLSWord(bn) & 1);
  481.     }
  482.  
  483.     /*
  484.      * This could be cached between calls to sieveBuild, but
  485.      * it's really not worth it; sieveSmall is *very* fast.
  486.      * sieveSmall returns a sieve of odd primes.
  487.      */
  488.     sieveSmall(small, SMALL);
  489.  
  490.     /*
  491.      * Okay, now sieve via the primes up to ssize*16+SMALLSTART-1,
  492.      * obtained from the small table.
  493.      */
  494.     i = (small[0] & 1) ? 0 : sieveSearch(small, SMALL, 0);
  495.     do {
  496.         p = 2 * i + SMALLSTART;
  497.  
  498.         /*
  499.          * Modulo is usually very expensive, but step is usually
  500.          * small, so this conditional is worth it.
  501.          */
  502.         t = (step < p) ? step : step % p;
  503.         if (!t) {
  504.             /*
  505.              * Instead of assert failing, returning all zero
  506.              * bits is the "correct" thing to do, but I think
  507.              * that the caller should take care of that
  508.              * themselves before starting.
  509.              */
  510.             assert(bnModQ(bn, p) != 0);
  511.             continue;
  512.         }
  513.         /*
  514.          * Get inverse of step mod p.  0 < t < p, and p is prime,
  515.          * so it has an inverse and sieveModInvert can't return 0.
  516.          */
  517.         t = sieveModInvert(t, p);
  518.         assert(t);
  519.         /* Negate t, so now t == -1/step (mod p) */
  520.         t = p - t;
  521.  
  522.         /* Now get the bignum modulo the prime. */
  523.         s = bnModQ(bn, p);
  524.  
  525.         /* Multiply by t, the negative inverse of step size */
  526. #if UINT_MAX/0xffff < 0xffff
  527.         s = (unsigned)(((unsigned long)s * t) % p);
  528. #else
  529.         s = (s * t) % p;
  530. #endif
  531.  
  532.         /* s is now the starting bit position, so sieve */
  533.         sieveSingle(array, size, s, p);
  534.  
  535.         /* Now do the double sieves as desired. */
  536.         for (j = 0; j < dbl; j++) {
  537.             /* Halve t modulo p */
  538. #if UINT_MAX < 0x1ffff
  539.             t = (t & 1) ? p/2 + t/2 + 1 : t/2;
  540.             /* Add t to s, modulo p with overflow checks. */
  541.             s += t;
  542.             if (s >= p || s < t)
  543.                 s -= p;
  544. #else
  545.             if (t & 1)
  546.                 t += p;
  547.             t /= 2;
  548.             /* Add t to s, modulo p */
  549.             s += t;
  550.             if (s >= p)
  551.                 s -= p;
  552. #endif
  553.             sieveSingle(array, size, s, p);
  554.         }
  555.  
  556.         /* And find the next prime */
  557.     } while ((i = sieveSearch(small, SMALL, i)) != 0);
  558.  
  559. #ifdef MSDOS
  560.     lbnMemFree(small, SMALL);
  561. #endif
  562.     return 0;    /* Success */
  563. }
  564.  
  565. /*
  566.  * Similar to the above, but use "step" (which must be even) as a step size
  567.  * rather than a fixed value of 2.  If "step" has any small divisors other
  568.  * than 2, this will blow up.
  569.  *
  570.  * Returns -1 on out of memory (MSDOS only, actually), and -2
  571.  * if step is found to be non-prime.
  572.  */
  573. int
  574. sieveBuildBig(unsigned char *array, unsigned size, struct BigNum const *bn,
  575.     struct BigNum const *step, unsigned dbl)
  576. {
  577.     unsigned i, j;    /* Loop index */
  578.     unsigned p;    /* Current small prime */
  579.     unsigned s;    /* Where to start operations in the big sieve */
  580.     unsigned t;    /* step modulo p, the current prime */
  581. #ifdef MSDOS    /* Use dynamic allocation rather than on the stack */
  582.     unsigned char *small;
  583. #else
  584.     unsigned char small[SMALL];
  585. #endif
  586.  
  587.     assert(array);
  588.  
  589. #ifdef MSDOS
  590.     small = lbnMemAlloc(SMALL);    /* Which allocator?  Not secure. */
  591.     if (!small)
  592.         return -1;    /* Failed */
  593. #endif
  594.     /*
  595.      * An odd step is a special case, since we must sieve by 2,
  596.      * which isn't in the small prime array and has a few other
  597.      * special properties.  These are:
  598.      * - Since the numbers are stored in binary, we don't need to
  599.      *   use bnModQ to find the remainder.
  600.      * - If step is odd, then t = step % 2 is 1, which allows
  601.      *   the elimination of a lot of math.  Inverting and negating
  602.      *   t don't change it, and multiplying s by 1 is a no-op,
  603.      *   so t isn't actually mentioned.
  604.      * - Since this is the first sieving, instead of calling
  605.      *   sieveSingle, we can just use memset to fill the array
  606.      *   with 0x55 or 0xAA.  Since a 1 bit means possible prime
  607.      *   (i.e. NOT divisible by 2), and the least significant bit
  608.      *   is first, if bn % 2 == 0, we use 0xAA (bit 0 = bn is NOT
  609.      *   prime), while if bn % 2 == 1, use 0x55.
  610.      *   (If step is even, bn must be odd, so fill the array with 0xFF.)
  611.      * - Any doublings need not be considered, since 2*bn+1 is odd, and
  612.      *   2*step is even, so none of these numbers are divisible by 2.
  613.      */
  614.     if (bnLSWord(step) & 1) {
  615.         s = bnLSWord(bn) & 1;
  616.         memset(array, 0xAA >> s, size);
  617.     } else {
  618.         /* Initialize the array to all 1's */
  619.         memset(array, 255, size);
  620.         assert(bnLSWord(bn) & 1);
  621.     }
  622.  
  623.     /*
  624.      * This could be cached between calls to sieveBuild, but
  625.      * it's really not worth it; sieveSmall is *very* fast.
  626.      * sieveSmall returns a sieve of the odd primes.
  627.      */
  628.     sieveSmall(small, SMALL);
  629.  
  630.     /*
  631.      * Okay, now sieve via the primes up to ssize*16+SMALLSTART-1,
  632.      * obtained from the small table.
  633.      */
  634.     i = (small[0] & 1) ? 0 : sieveSearch(small, SMALL, 0);
  635.     do {
  636.         p = 2 * i + SMALLSTART;
  637.  
  638.         t = bnModQ(step, p);
  639.         if (!t) {
  640.             assert(bnModQ(bn, p) != 0);
  641.             continue;
  642.         }
  643.         /* Get negative inverse of step */
  644.         t = sieveModInvert(bnModQ(step, p), p);
  645.         assert(t);
  646.         t = p-t;
  647.  
  648.         /* Okay, we have a prime - get the remainder */
  649.         s = bnModQ(bn, p);
  650.  
  651.         /* Now multiply s by the negative inverse of step (mod p) */
  652. #if UINT_MAX < 0xffff * 0xffff
  653.         s = (unsigned)(((unsigned long)s * t) % p);
  654. #else
  655.         s = (s * t) % p;
  656. #endif
  657.         /* We now have the starting bit pos */
  658.         sieveSingle(array, size, s, p);
  659.  
  660.         /* Now do the double sieves as desired. */
  661.         for (j = 0; j < dbl; j++) {
  662.             /* Halve t modulo p */
  663. #if UINT_MAX < 0x1ffff
  664.             t = (t & 1) ? p/2 + t/2 + 1 : t/2;
  665.             /* Add t to s, modulo p with overflow checks. */
  666.             s += t;
  667.             if (s >= p || s < t)
  668.                 s -= p;
  669. #else
  670.             if (t & 1)
  671.                 t += p;
  672.             t /= 2;
  673.             /* Add t to s, modulo p */
  674.             s += t;
  675.             if (s >= p)
  676.                 s -= p;
  677. #endif
  678.             sieveSingle(array, size, s, p);
  679.         }
  680.  
  681.         /* And find the next prime */
  682.     } while ((i = sieveSearch(small, SMALL, i)) != 0);
  683.  
  684. #ifdef MSDOS
  685.     lbnMemFree(small, SMALL);
  686. #endif
  687.     return 0;    /* Success */
  688. }
  689.