home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / crypl200.zip / BNLIB / PRIME.C < prev    next >
C/C++ Source or Header  |  1996-05-16  |  21KB  |  650 lines
  1. /*
  2.  * Prime generation using the bignum library and sieving.
  3.  *
  4.  * Copyright (c) 1995  Colin Plumb.  All rights reserved.
  5.  * For licensing and other legal details, see the file legal.c.
  6.  */
  7. #ifndef HAVE_CONFIG_H
  8. #define HAVE_CONFIG_H 0
  9. #endif
  10. #if HAVE_CONFIG_H
  11. #include "config.h"
  12. #endif
  13.  
  14. /*
  15.  * Some compilers complain about #if FOO if FOO isn't defined,
  16.  * so do the ANSI-mandated thing explicitly...
  17.  */
  18. #ifndef NO_ASSERT_H
  19. #define NO_ASSERT_H 0
  20. #endif
  21. #if !NO_ASSERT_H
  22. #include <assert.h>
  23. #else
  24. #define assert(x) (void)0
  25. #endif
  26.  
  27. #include <stdarg.h>    /* We just can't live without this... */
  28.  
  29. #ifndef BNDEBUG
  30. #define BNDEBUG 0
  31. #endif
  32. #if BNDEBUG
  33. #include <stdio.h>
  34. #endif
  35.  
  36. #include "bn.h"
  37. #include "lbnmem.h"
  38. #include "prime.h"
  39. #include "sieve.h"
  40.  
  41. #include "kludge.h"
  42.  
  43. /* Size of the shuffle table */
  44. #define SHUFFLE    256
  45. /* Size of the sieve area */
  46. #define SIEVE 32768u/16
  47.  
  48. /*
  49.  * Helper function that does the slow primality test.
  50.  * bn is the input bignum; a and e are temporary buffers that are
  51.  * allocated by the caller to save overhead.
  52.  *
  53.  * Returns 0 if prime, >0 if not prime, and -1 on error (out of memory).
  54.  * If not prime, returns the number of modular exponentiations performed.
  55.  * Calls the fiven progress function with a '*' for each primality test
  56.  * that is passed.
  57.  *
  58.  * The testing consists of strong pseudoprimality tests, to the bases
  59.  * 2, 3, 5, 7, 11, 13 and 17.  (Also called Miller-Rabin, although that's
  60.  * not technically correct if we're using fixed bases.)  Some people worry
  61.  * that this might not be enough.  Number theorists may wish to generate
  62.  * primality proofs, but for random inputs, this returns non-primes with a
  63.  * probability which is quite negligible, which is good enough.
  64.  *
  65.  * It has been proved (see Carl Pomerance, "On the Distribution of
  66.  * Pseudoprimes", Math. Comp. v.37 (1981) pp. 587-593) that the number
  67.  * of pseudoprimes (composite numbers that pass a Fermat test to the
  68.  * base 2) less than x is bounded by:
  69.  * exp(ln(x)^(5/14)) <= P_2(x)    ### CHECK THIS FORMULA - it looks wrong! ###
  70.  * P_2(x) <= x * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))).
  71.  * Thus, the local density of Pseudoprimes near x is at most
  72.  * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))), and at least
  73.  * exp(ln(x)^(5/14) - ln(x)).  Here are some values of this function
  74.  * for various k-bit numbers x = 2^k:
  75.  * Bits    Density <=    Bit equivalent    Density >=    Bit equivalent
  76.  *  128    3.577869e-07     21.414396    4.202213e-37     120.840190
  77.  *  192    4.175629e-10     31.157288    4.936250e-56     183.724558
  78.  *  256 5.804314e-13     40.647940    4.977813e-75     246.829095
  79.  *  384 1.578039e-18     59.136573    3.938861e-113     373.400096
  80.  *  512 5.858255e-24     77.175803    2.563353e-151     500.253110
  81.  *  768 1.489276e-34    112.370944    7.872825e-228     754.422724
  82.  * 1024 6.633188e-45    146.757062    1.882404e-304    1008.953565
  83.  *
  84.  * As you can see, there's quite a bit of slop between these estimates.
  85.  * In fact, the density of pseudoprimes is conjectured to be closer
  86.  * to the square of that upper bound.  E.g. the density of pseudoprimes
  87.  * of size 256 is around 3 * 10^-27.  The density of primes is very
  88.  * high, from 0.005636 at 256 bits to 0.001409 at 1024 bits, i.e.
  89.  * more than 10^-3.
  90.  *
  91.  * For those people used to cryptographic levels of security where the
  92.  * 56 bits of DES key space is too small because it's exhaustible with
  93.  * custom hardware searching engines, note that you are not generating
  94.  * 50,000,000 primes per second on each of 56,000 custom hardware chips
  95.  * for several hours.  The chances that another Dinosaur Killer asteroid
  96.  * will land today is about 10^-11 or 2^-36, so it would be better to
  97.  * spend your time worrying about *that*.  Well, okay, there should be
  98.  * some derating for the chance that astronomers haven't seen it yet, but
  99.  * I think you get the idea.  For a good feel about the probability of
  100.  * various events, I have heard that a good book is by E'mile Borel,
  101.  * "Les Probabilite's et la vie".  (The 's are accents, not apostrophes.)
  102.  *
  103.  * For more on the subject, try "Finding Four Million Large Random Primes",
  104.  * by Ronald Rivest, in Advancess in Cryptology: Proceedings of Crypto '90.
  105.  * He used a small-divisor test, then a Fermat test to the base 2, and then
  106.  * 8 iterations of a Miller-Rabin test.  About 718 million random 256-bit
  107.  * integers were generated, 43,741,404 passed the small divisor test,
  108.  * 4,058,000 passed the Fermat test, and all 4,058,000 passed all 8
  109.  * iterations of the Miller-Rabin test, proving their primality beyond most
  110.  * reasonable doubts.
  111.  *
  112.  * If the probability of getting a pseudoprime is some small p, then
  113.  * the probability of not getting it in t trials is (1-p)^t.  Remember
  114.  * that, for small p, (1-p)^(1/p) ~ 1/e, the base of natural logarithms.
  115.  * (This is more commonly expressed as e = lim_{x\to\infty} (1+1/x)^x.)
  116.  * Thus, (1-p)^t ~ e^(-p*t) = exp(-p*t).  So the odds of being able to
  117.  * do this many tests without seeing a pseudoprime if you assume that
  118.  * p = 10^-6 (one in a million) is one in 57.86.  If you assume that
  119.  * p = 2*10^-6, it's one in 3347.6.  So it's implausible that the
  120.  * density of pseudoprimes is much more than one millionth the density
  121.  * of primes.
  122.  *
  123.  * He also gives a theoretical argument that the chance of finding a
  124.  * 256-bit non-prime which satisfies one Fermat test to the base 2 is less
  125.  * than 10^-22.  The small divisor test improves this number, and if the
  126.  * numbers are 512 bits (as needed for a 1024-bit key) the odds of failure
  127.  * shrink to about 10^-44.  Thus, he concludes, for practical purposes
  128.  * *one* Fermat test to the base 2 is sufficient.
  129.  */
  130. static int
  131. primeTest(struct BigNum const *bn, struct BigNum *e, struct BigNum *a,
  132.     int (*f)(void *arg, int c), void *arg)
  133. {
  134.     unsigned i, j;
  135.     unsigned k, l;
  136.     int err;
  137.     static unsigned const primes[] = {2, 3, 5, 7, 11, 13, 17};
  138.  
  139. #if BNDEBUG    /* Debugging */
  140.     /*
  141.      * This is debugging code to test the sieving stage.
  142.      * If the sieving is wrong, it will let past numbers with
  143.      * small divisors.  The prime test here will still work, and
  144.      * weed them out, but you'll be doing a lot more slow tests,
  145.      * and presumably excluding from consideration some other numbers
  146.      * which might be prime.  This check just verifies that none
  147.      * of the candidates have any small divisors.  If this
  148.      * code is enabled and never triggers, you can feel quite
  149.      * confident that the sieving is doing its job.
  150.      */
  151.     i = bnModQ(bn, 30030);    /* 30030 = 2 * 3 * 5 * 7 * 11 * 13 */
  152.     if (!(i % 2)) printf("bn div by 2!");
  153.     if (!(i % 3)) printf("bn div by 3!");
  154.     if (!(i % 5)) printf("bn div by 5!");
  155.     if (!(i % 7)) printf("bn div by 7!");
  156.     if (!(i % 11)) printf("bn div by 11!");
  157.     if (!(i % 13)) printf("bn div by 13!");
  158.     i = bnModQ(bn, 7429);    /* 7429 = 17 * 19 * 23 */
  159.     if (!(i % 17)) printf("bn div by 17!");
  160.     if (!(i % 19)) printf("bn div by 19!");
  161.     if (!(i % 23)) printf("bn div by 23!");
  162.     i = bnModQ(bn, 33263);    /* 33263 = 29 * 31 * 37 */
  163.     if (!(i % 29)) printf("bn div by 29!");
  164.     if (!(i % 31)) printf("bn div by 31!");
  165.     if (!(i % 37)) printf("bn div by 37!");
  166. #endif
  167.  
  168.     /*
  169.      * Now, check that bn is prime.  If it passes to the base 2,
  170.      * it's prime beyond all reasonable doubt, and everything else
  171.      * is just gravy, but it gives people warm fuzzies to do it.
  172.      *
  173.      * This starts with verifying Euler's criterion for a base of 2.
  174.      * This is the fastest pseudoprimality test that I know of,
  175.      * saving a modular squaring over a Fermat test, as well as
  176.      * being stronger.  7/8 of the time, it's as strong as a strong
  177.      * pseudoprimality test, too.  (The exception being when bn ==
  178.      * 1 mod 8 and 2 is a quartic residue, i.e. bn is of the form
  179.      * a^2 + (8*b)^2.)  The precise series of tricks used here is
  180.      * not documented anywhere, so here's an explanation.
  181.      * Euler's criterion states that if p is prime then a^((p-1)/2)
  182.      * is congruent to Jacobi(a,p), modulo p.  Jacobi(a,p) is
  183.      * a function which is +1 if a is a square modulo p, and -1 if
  184.      * it is not.  For a = 2, this is particularly simple.  It's
  185.      * +1 if p == +/-1 (mod 8), and -1 if m == +/-3 (mod 8).
  186.      * If p == 3 mod 4, then all a strong test does is compute
  187.      * 2^((p-1)/2). and see if it's +1 or -1.  (Euler's criterion
  188.      * says *which* it should be.)  If p == 5 (mod 8), then
  189.      * 2^((p-1)/2) is -1, so the initial step in a strong test,
  190.      * looking at 2^((p-1)/4), is wasted - you're not going to
  191.      * find a +/-1 before then if it *is* prime, and it shouldn't
  192.      * have either of those values if it isn't.  So don't bother.
  193.      *
  194.      * The remaining case is p == 1 (mod 8).  In this case, we
  195.      * expect 2^((p-1)/2) == 1 (mod p), so we expect that the
  196.      * square root of this, 2^((p-1)/4), will be +/-1 (mod p).
  197.      * Evaluating this saves us a modular squaring 1/4 of the time.
  198.      * If it's -1, a strong pseudoprimality test would call p
  199.      * prime as well.  Only if the result is +1, indicating that
  200.      * 2 is not only a quadratic residue, but a quartic one as well,
  201.      * does a strong pseudoprimality test verify more things than
  202.      * this test does.  Good enough.
  203.      *
  204.      * We could back that down another step, looking at 2^((p-1)/8)
  205.      * if there was a cheap way to determine if 2 were expected to
  206.      * be a quartic residue or not.  Dirichlet proved that 2 is
  207.      * a quartic residue iff p is of the form a^2 + (8*b^2).
  208.      * All primes == 1 (mod 4) can be expressed as a^2 + (2*b)^2,
  209.      * but I see no cheap way to evaluate this condition.
  210.      */
  211.     if (bnCopy(e, bn) < 0)
  212.         return -1;
  213.     (void)bnSubQ(e, 1);
  214.     l = bnLSWord(e);
  215.  
  216.     j = 1;    /* Where to start in prime array for strong prime tests */
  217.  
  218.     if (l & 7) {
  219.         bnRShift(e, 1);
  220.         if (bnTwoExpMod(a, e, bn) < 0)
  221.             return -1;
  222.         if ((l & 7) == 6) {
  223.             /* bn == 7 mod 8, expect +1 */
  224.             if (bnBits(a) != 1)
  225.                 return 1;    /* Not prime */
  226.             k = 1;
  227.         } else {
  228.             /* bn == 3 or 5 mod 8, expect -1 == bn-1 */
  229.             if (bnAddQ(a, 1) < 0)
  230.                 return -1;
  231.             if (bnCmp(a, bn) != 0)
  232.                 return 1;    /* Not prime */
  233.             k = 1;
  234.             if (l & 4) {
  235.                 /* bn == 5 mod 8, make odd for strong tests */
  236.                 bnRShift(e, 1);
  237.                 k = 2;
  238.             }
  239.         }
  240.     } else {
  241.         /* bn == 1 mod 8, expect 2^((bn-1)/4) == +/-1 mod bn */
  242.         bnRShift(e, 2);
  243.         if (bnTwoExpMod(a, e, bn) < 0)
  244.             return -1;
  245.         if (bnBits(a) == 1) {
  246.             j = 0;    /* Re-do strong prime test to base 2 */
  247.         } else {
  248.             if (bnAddQ(a, 1) < 0)
  249.                 return -1;
  250.             if (bnCmp(a, bn) != 0)
  251.                 return 1;    /* Not prime */
  252.         }
  253.         k = 2 + bnMakeOdd(e);
  254.     }
  255.     /* It's prime!  Now go on to confirmation tests */
  256.  
  257.     /*
  258.      * Now, e = (bn-1)/2^k is odd.  k >= 1, and has a given value
  259.      * with probability 2^-k, so its expected value is 2.
  260.      * j = 1 in the usual case when the previous test was as good as
  261.      * a strong prime test, but 1/8 of the time, j = 0 because
  262.      * the strong prime test to the base 2 needs to be re-done.
  263.      */
  264.     for (i = j; i < sizeof(primes)/sizeof(*primes); i++) {
  265.         if (f && (err = f(arg, '*')) < 0)
  266.             return err;
  267.         (void)bnSetQ(a, primes[i]);
  268.         if (bnExpMod(a, a, e, bn) < 0)
  269.             return -1;
  270.         if (bnBits(a) == 1)
  271.             continue;    /* Passed this test */
  272.  
  273.         l = k;
  274.         for (;;) {
  275.             if (bnAddQ(a, 1) < 0)
  276.                 return -1;
  277.             if (bnCmp(a, bn) == 0)    /* Was result bn-1? */
  278.                 break;    /* Prime */
  279.             if (!--l)    /* Reached end, not -1? luck? */
  280.                 return i+2-j;    /* Failed, not prime */
  281.             /* This portion is executed, on average, once. */
  282.             (void)bnSubQ(a, 1);    /* Put a back where it was. */
  283.             if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0)
  284.                 return -1;
  285.             if (bnBits(a) == 1)
  286.                 return i+2-j;    /* Failed, not prime */
  287.         }
  288.         /* It worked (to the base primes[i]) */
  289.     }
  290.     
  291.     /* Yes, we've decided that it's prime. */
  292.     if (f && (err = f(arg, '*')) < 0)
  293.         return err;
  294.     return 0;    /* Prime! */
  295. }
  296.  
  297. /*
  298.  * Modifies the bignum to return a nearby (slightly larger) number which
  299.  * is a probable prime.  Returns >=0 on success or -1 on failure (out of
  300.  * memory).  The return value is the number of unsuccessful modular
  301.  * exponentiations performed.  This never gives up searching.
  302.  *
  303.  * All other arguments are optional.  They may be NULL.  They are:
  304.  *
  305.  * unsigned (*rand)(unsigned limit)
  306.  * For better distributed numbers, supply a non-null pointer to a
  307.  * function which returns a random x, 0 <= x < limit.  (It may make it
  308.  * simpler to know that 0 < limit <= SHUFFLE, so you need at most a byte.)
  309.  * The program generates a large window of sieve data and then does
  310.  * pseudoprimality tests on the data.  If a rand function is supplied,
  311.  * the candidates which survive sieving are shuffled with a window of
  312.  * size SHUFFLE before testing to increase the uniformity of the prime
  313.  * selection.  This isn't perfect, but it reduces the correlation between
  314.  * the size of the prime-free gap before a prime and the probability
  315.  * that that prime will be found by a sequential search.
  316.  *
  317.  * If rand is NULL, sequential search is used.  If you want sequential
  318.  * search, note that the search begins with the given number; if you're
  319.  * trying to generate consecutive primes, you must increment the previous
  320.  * one by two before calling this again.
  321.  *
  322.  * int (*f)(void *arg, int c), void *arg
  323.  * The function f argument, if non-NULL, is called with progress indicator
  324.  * characters for printing.  A dot (.) is written every time a primality test
  325.  * is failed, a star (*) every time one is passed, and a slash (/) in the
  326.  * (very rare) case that the sieve was emptied without finding a prime
  327.  * and is being refilled.  f is also passed the void *arg argument for
  328.  * private context storage.  If f returns < 0, the test aborts and returns
  329.  * that value immediately.
  330.  *
  331.  * The "exponent" argument, and following unsigned numbers, are exponents
  332.  * for which an inverse is desired, modulo p.  For a d to exist such that
  333.  * (x^e)^d == x (mod p), then d*e == 1 (mod p-1), so gcd(e,p-1) must be 1.
  334.  * The prime returned is constrained to not be congruent to 1 modulo
  335.  * any of the zero-terminated list of 16-bit numbers.  Note that this list
  336.  * should contain all the small prime factors of e.  (You'll have to test
  337.  * for large prime factors of e elsewhere, but the chances of needing to
  338.  * generate another prime are low.)
  339.  *
  340.  * The list is terminated by a 0, and may be empty.
  341.  *
  342.  */
  343. int
  344. primeGen(struct BigNum *bn, unsigned (*rand)(unsigned),
  345.          int (*f)(void *arg, int c), void *arg, unsigned exponent, ...)
  346. {
  347.     int retval;
  348.     int modexps = 0;
  349.     unsigned short offsets[SHUFFLE];
  350.     unsigned i, j;
  351.     unsigned p, q, prev;
  352.     struct BigNum a, e;
  353. #ifdef MSDOS
  354.     unsigned char *sieve;
  355. #else
  356.     unsigned char sieve[SIEVE];
  357. #endif
  358.  
  359. #ifdef MSDOS
  360.     sieve = lbnMemAlloc(SIEVE);
  361.     if (!sieve)
  362.         return -1;
  363. #endif
  364.  
  365.     bnBegin(&a);
  366.     bnBegin(&e);
  367.  
  368. #if 0    /* Self-test (not used for production) */
  369. {
  370.     struct BigNum t;
  371.     static unsigned char const prime1[] = {5};
  372.     static unsigned char const prime2[] = {7};
  373.     static unsigned char const prime3[] = {11};
  374.     static unsigned char const prime4[] = {1, 1}; /* 257 */
  375.     static unsigned char const prime5[] = {0xFF, 0xF1}; /* 65521 */
  376.     static unsigned char const prime6[] = {1, 0, 1}; /* 65537 */
  377.     static unsigned char const prime7[] = {1, 0, 3}; /* 65539 */
  378.     /* A small prime: 1234567891 */
  379.     static unsigned char const prime8[] = {0x49, 0x96, 0x02, 0xD3};
  380.     /* A slightly larger prime: 12345678901234567891 */
  381.     static unsigned char const prime9[] = {
  382.         0xAB, 0x54, 0xA9, 0x8C, 0xEB, 0x1F, 0x0A, 0xD3 };
  383.     /*
  384.      * No, 123456789012345678901234567891 isn't prime; it's just a
  385.      * lucky, easy-to-remember conicidence.  (You have to go to
  386.      * ...4567907 for a prime.)
  387.      */
  388.     static struct {
  389.         unsigned char const *prime;
  390.         unsigned size;
  391.     } const primelist[] = {
  392.         { prime1, sizeof(prime1) },
  393.         { prime2, sizeof(prime2) },
  394.         { prime3, sizeof(prime3) },
  395.         { prime4, sizeof(prime4) },
  396.         { prime5, sizeof(prime5) },
  397.         { prime6, sizeof(prime6) },
  398.         { prime7, sizeof(prime7) },
  399.         { prime8, sizeof(prime8) },
  400.         { prime9, sizeof(prime9) } };
  401.  
  402.     bnBegin(&t);
  403.  
  404.     for (i = 0; i < sizeof(primelist)/sizeof(primelist[0]); i++) {
  405.             bnInsertBytes(&t, primelist[i].prime, 0,
  406.                       primelist[i].size);
  407.             bnCopy(&e, &t);
  408.             (void)bnSubQ(&e, 1);
  409.             bnTwoExpMod(&a, &e, &t);
  410.             p = bnBits(&a);
  411.             if (p != 1) {
  412.                 printf(
  413.             "Bug: Fermat(2) %u-bit output (1 expected)\n", p);
  414.                 fputs("Prime = 0x", stdout);
  415.                 for (j = 0; j < primelist[i].size; j++)
  416.                     printf("%02X", primelist[i].prime[j]);
  417.                 putchar('\n');
  418.             }
  419.             bnSetQ(&a, 3);
  420.             bnExpMod(&a, &a, &e, &t);
  421.             if (p != 1) {
  422.                 printf(
  423.             "Bug: Fermat(3) %u-bit output (1 expected)\n", p);
  424.                 fputs("Prime = 0x", stdout);
  425.                 for (j = 0; j < primelist[i].size; j++)
  426.                     printf("%02X", primelist[i].prime[j]);
  427.                 putchar('\n');
  428.             }
  429.         }
  430.  
  431.     bnEnd(&t);
  432. }
  433. #endif
  434.  
  435.     /* First, make sure that bn is odd. */
  436.     if ((bnLSWord(bn) & 1) == 0)
  437.         (void)bnAddQ(bn, 1);
  438.  
  439. retry:
  440.     /* Then build a sieve starting at bn. */
  441.     sieveBuild(sieve, SIEVE, bn, 2, 0);
  442.  
  443.     /* Do the extra exponent sieving */
  444.     if (exponent) {
  445.         va_list ap;
  446.         unsigned t = exponent;
  447.  
  448.         va_start(ap, exponent);
  449.  
  450.         do {
  451.             /* The exponent had better be odd! */
  452.             assert(t & 1);
  453.  
  454.             i = bnModQ(bn, t);
  455.             /* Find 1-i */
  456.             if (i == 0)
  457.                 i = 1;
  458.             else if (--i)
  459.                 i = t - i;
  460.  
  461.             /* Divide by 2, modulo the exponent */
  462.             i = (i & 1) ? i/2 + t/2 + 1 : i/2;
  463.  
  464.             /* Remove all following multiples from the sieve. */
  465.             sieveSingle(sieve, SIEVE, i, t);
  466.  
  467.             /* Get the next exponent value */
  468.             t = va_arg(ap, unsigned);
  469.         } while (t);
  470.  
  471.         va_end(ap);
  472.     }
  473.  
  474.     /* Fill up the offsets array with the first SHUFFLE candidates */
  475.     i = p = 0;
  476.     /* Get first prime */
  477.     if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0)
  478.         offsets[i++] = p;
  479.     /*
  480.      * If we have a prime and we have a shuffle function so it makes
  481.      * to fill up the table, so do.
  482.      */
  483.     if (rand && i) {
  484.         do {
  485.             p = sieveSearch(sieve, SIEVE, p);
  486.             if (p == 0)
  487.                 break;
  488.             offsets[i++] = p;
  489.         } while (i < SHUFFLE);
  490.     }
  491.  
  492.     /* Choose a random candidate for experimentation */
  493.     prev = 0;
  494.     while (i) {
  495.         /* Pick a random entry from the shuffle table */
  496.         j = rand ? rand(i) : 0;
  497.         q = offsets[j];
  498.  
  499.         /* Replace the entry with some more data, if possible */
  500.         if (p && (p = sieveSearch(sieve, SIEVE, p)) != 0)
  501.             offsets[j] = p;
  502.         else
  503.             offsets[j] = offsets[--i];
  504.  
  505.         /* Adjust bn to have the right value */
  506.         if (q > prev) {
  507.             if (bnAddQ(bn, q-prev) < 0)
  508.                 goto failed;
  509.             if (bnAddQ(bn, q-prev) < 0)
  510.                 goto failed;
  511.         } else {
  512.             if (bnSubQ(bn, prev-q))
  513.                 goto failed;
  514.             if (bnSubQ(bn, prev-q))
  515.                 goto failed;
  516.         }
  517.         prev = q;
  518.  
  519.         /* Now do the Fermat tests */
  520.         retval = primeTest(bn, &e, &a, f, arg);
  521.         if (retval <= 0)
  522.             goto done;    /* Success or error */
  523.         modexps += retval;
  524.         if (f && (retval = f(arg, '.')) < 0)
  525.             goto done;
  526.     }
  527.  
  528.     /* Ran out of sieve space - increase bn and keep trying. */
  529.  
  530.     if (bnSubQ(bn, (SIEVE*8)-prev))
  531.         goto failed;
  532.     if (bnSubQ(bn, (SIEVE*8)-prev))
  533.         goto failed;
  534.     if (f && (retval = f(arg, '/')) < 0)
  535.         goto done;
  536.     goto retry;
  537.  
  538. failed:
  539.     retval = -1;
  540. done:
  541.     bnEnd(&e);
  542.     bnEnd(&a);
  543.     lbnMemWipe(offsets, sizeof(offsets));
  544. #ifdef MSDOS
  545.     lbnMemFree(sieve, SIEVE);
  546. #else
  547.     lbnMemWipe(sieve, sizeof(sieve));
  548. #endif
  549.  
  550.     return retval < 0 ? retval : modexps;
  551. }
  552.  
  553. /*
  554.  * Similar, but searches forward from the given starting value in steps of
  555.  * "step" rather than 1.  The step size must be even, and bn must be odd.
  556.  * Among other possibilities, this can be used to generate "strong"
  557.  * primes, where p-1 has a large prime factor.
  558.  */
  559. int
  560. primeGenStrong(struct BigNum *bn, struct BigNum const *step,
  561.     int (*f)(void *arg, int c), void *arg)
  562. {
  563.     int retval;
  564.     unsigned p, prev;
  565.     struct BigNum a, e;
  566.     int modexps = 0;
  567. #ifdef MSDOS
  568.     unsigned char *sieve;
  569. #else
  570.     unsigned char sieve[SIEVE];
  571. #endif
  572.  
  573. #ifdef MSDOS
  574.     sieve = lbnMemAlloc(SIEVE);
  575.     if (!sieve)
  576.         return -1;
  577. #endif
  578.  
  579.     /* Step must be even and bn must be odd */
  580.     assert((bnLSWord(step) & 1) == 0);
  581.     assert((bnLSWord(bn) & 1) == 1);
  582.  
  583.     bnBegin(&a);
  584.     bnBegin(&e);
  585.  
  586.     for (;;) {
  587.         if (sieveBuildBig(sieve, SIEVE, bn, step, 0) < 0)
  588.             goto failed;
  589.  
  590.         p = prev = 0;
  591.         if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) {
  592.             do {
  593.                 /*
  594.                  * Adjust bn to have the right value,
  595.                  * adding (p-prev) * 2*step.
  596.                  */
  597.                 assert(p >= prev);
  598.                 /* Compute delta into a */
  599.                 if (bnMulQ(&a, step, p-prev) < 0)
  600.                     goto failed;
  601.                 if (bnAdd(bn, &a) < 0)
  602.                     goto failed;
  603.                 prev = p;
  604.  
  605.                 retval = primeTest(bn, &e, &a, f, arg);
  606.                 if (retval <= 0)
  607.                     goto done;    /* Success! */
  608.                 modexps += retval;
  609.                 if (f && (retval = f(arg, '.')) < 0)
  610.                     goto done;
  611.  
  612.                 /* And try again */
  613.                 p = sieveSearch(sieve, SIEVE, p);
  614.             } while (p);
  615.         }
  616.  
  617.         /* Ran out of sieve space - increase bn and keep trying. */
  618. #if SIEVE*8 == 65536
  619.         /* Corner case that will never actually happen */
  620.         if (!prev) {
  621.             if (bnAdd(bn, step) < 0)
  622.                 goto failed;
  623.             p = 65535;
  624.         } else {
  625.             p = (unsigned)(SIEVE*8 - prev);
  626.         }
  627. #else
  628.         p = SIEVE*8 - prev;
  629. #endif
  630.         if (bnMulQ(&a, step, p) < 0 || bnAdd(bn, &a) < 0)
  631.             goto failed;
  632.         if (f && (retval = f(arg, '/')) < 0)
  633.             goto done;
  634.     } /* for (;;) */
  635.  
  636. failed:
  637.     retval = -1;
  638.  
  639. done:
  640.  
  641.     bnEnd(&e);
  642.     bnEnd(&a);
  643. #ifdef MSDOS
  644.     lbnMemFree(sieve, SIEVE);
  645. #else
  646.     lbnMemWipe(sieve, sizeof(sieve));
  647. #endif
  648.     return retval < 0 ? retval : modexps;
  649. }
  650.