home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / crypl200.zip / BNLIB / GERMAIN.C < prev    next >
C/C++ Source or Header  |  1996-05-16  |  16KB  |  469 lines
  1. /*
  2.  * Sophie Germain 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. #ifndef BNDEBUG
  28. #define BNDEBUG 0
  29. #endif
  30. #if BNDEBUG
  31. #include <stdio.h>
  32. #endif
  33.  
  34. #include "bn.h"
  35. #include "germain.h"
  36. #include "jacobi.h"
  37. #include "lbnmem.h"    /* For lbnMemWipe */
  38. #include "sieve.h"
  39.  
  40. #include "kludge.h"
  41.  
  42. /* Size of the sieve area (can be up to 65536/8 = 8192) */
  43. #define SIEVE 8192
  44.  
  45. /*
  46.  * Helper function that does the slow primality test.
  47.  * bn is the input bignum; a and e are temporary buffers that are
  48.  * allocated by the caller to save overhead.  bn2 is filled with
  49.  * a copy of 2*bn+1 if bn is found to be prime.
  50.  *
  51.  * Returns 0 if both bn abd bn2 are prime, >0 if not prime, and -1 on
  52.  * error (out of memory).  If not prime, the return value is the number
  53.  * of modular exponentiations performed.   Prints a '+' or '-' on the
  54.  * given FILE (if any) for each test that is passed by bn, and a '*'
  55.  * for each test that is passed by bn2.
  56.  *
  57.  * The testing consists of strong pseudoprimality tests, to the bases 2,
  58.  * 3, 5, 7, 11, 13 and 17.  (Also called Miller-Rabin, although that's not
  59.  * technically correct if we're using fixed bases.)  Some people worry
  60.  * that this might not be enough.  Number theorists may wish to generate
  61.  * primality proofs, but for random inputs, this returns non-primes with a
  62.  * probability which is quite negligible, which is good enough.
  63.  *
  64.  * It has been proved (see Carl Pomerance, "On the Distribution of
  65.  * Pseudoprimes", Math. Comp. v.37 (1981) pp. 587-593) that the number
  66.  * of pseudoprimes (composite numbers that pass a Fermat test to the
  67.  * base 2) less than x is bounded by:
  68.  * exp(ln(x)^(5/14)) <= P_2(x)    ### CHECK THIS FORMULA - it looks wrong! ###
  69.  * P_2(x) <= x * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))).
  70.  * Thus, the local density of Pseudoprimes near x is at most
  71.  * exp(-1/2 * ln(x) * ln(ln(ln(x))) / ln(ln(x))), and at least
  72.  * exp(ln(x)^(5/14) - ln(x)).  Here are some values of this function
  73.  * for various k-bit numbers x = 2^k:
  74.  * Bits    Density <=    Bit equivalent    Density >=    Bit equivalent
  75.  *  128    3.577869e-07     21.414396    4.202213e-37     120.840190
  76.  *  192    4.175629e-10     31.157288    4.936250e-56     183.724558
  77.  *  256 5.804314e-13     40.647940    4.977813e-75     246.829095
  78.  *  384 1.578039e-18     59.136573    3.938861e-113     373.400096
  79.  *  512 5.858255e-24     77.175803    2.563353e-151     500.253110
  80.  *  768 1.489276e-34    112.370944    7.872825e-228     754.422724
  81.  * 1024 6.633188e-45    146.757062    1.882404e-304    1008.953565
  82.  *
  83.  * As you can see, there's quite a bit of slop between these estimates.
  84.  * In fact, the density of pseudoprimes is conjectured to be closer
  85.  * to the square of that upper bound.  E.g. the density of pseudoprimes
  86.  * of size 256 is around 3 * 10^-27.  The density of primes is very
  87.  * high, from 0.005636 at 256 bits to 0.001409 at 1024 bits, i.e.
  88.  * more than 10^-3.
  89.  *
  90.  * For those people used to cryptographic levels of security where the
  91.  * 56 bits of DES key space is too small because it's exhaustible with
  92.  * custom hardware searching engines, note that you are not generating
  93.  * 50,000,000 primes per second on each of 56,000 custom hardware chips
  94.  * for several hours.  The chances that another Dinosaur Killer asteroid
  95.  * will land today is about 10^-11 or 2^-36, so it would be better to
  96.  * spend your time worrying about *that*.  Well, okay, there should be
  97.  * some derating for the chance that astronomers haven't seen it yet, but
  98.  * I think you get the idea.  For a good feel about the probability of
  99.  * various events, I have heard that a good book is by E'mile Borel,
  100.  * "Les Probabilite's et la vie".  (The 's are accents, not apostrophes.)
  101.  *
  102.  * For more on the subject, try "Finding Four Million Large Random Primes",
  103.  * by Ronald Rivest, in Advancess in Cryptology: Proceedings of Crypto '90.
  104.  * He used a small-divisor test, then a Fermat test to the base 2, and then
  105.  * 8 iterations of a Miller-Rabin test.  About 718 million random 256-bit
  106.  * integers were generated, 43,741,404 passed the small divisor test,
  107.  * 4,058,000 passed the Fermat test, and all 4,058,000 passed all 8
  108.  * iterations of the Miller-Rabin test, proving their primality beyond most
  109.  * reasonable doubts.
  110.  *
  111.  * If the probability of getting a pseudoprime is some small p, then
  112.  * the probability of not getting it in t trials is (1-p)^t.  Remember
  113.  * that, for small p, (1-p)^(1/p) ~ 1/e, the base of natural logarithms.
  114.  * (This is more commonly expressed as e = lim_{x\to\infty} (1+1/x)^x.)
  115.  * Thus, (1-p)^t ~ e^(-p*t) = exp(-p*t).  So the odds of being able to
  116.  * do this many tests without seeing a pseudoprime if you assume that
  117.  * p = 10^-6 (one in a million) is one in 57.86.  If you assume that
  118.  * p = 2*10^-6, it's one in 3347.6.  So it's implausible that the
  119.  * density of pseudoprimes is much more than one millionth the density
  120.  * of primes.
  121.  *
  122.  * He also gives a theoretical argument that the chance of finding a
  123.  * 256-bit non-prime which satisfies one Fermat test to the base 2 is less
  124.  * than 10^-22.  The small divisor test improves this number, and if the
  125.  * numbers are 512 bits (as needed for a 1024-bit key) the odds of failure
  126.  * shrink to about 10^-44.  Thus, he concludes, for practical purposes
  127.  * *one* Fermat test to the base 2 is sufficient.
  128.  */
  129. static int
  130. germainPrimeTest(struct BigNum const *bn, struct BigNum *bn2, struct BigNum *e,
  131.     struct BigNum *a, int (*f)(void *arg, int c), void *arg)
  132. {
  133.     int err;
  134.     unsigned i;
  135.     int j;
  136.     unsigned k, l;
  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, 15015);    /* 15015 = 3 * 5 * 7 * 11 * 13 */
  152.     if (!(i % 3)) printf("bn div by 3!");
  153.     if ((i % 3) == 1) printf("bn2 div by 3!");
  154.     if (!(i % 5)) printf("bn div by 5!");
  155.     if ((i % 5) == 2) printf("bn2 div by 5!");
  156.     if (!(i % 7)) printf("bn div by 7!");
  157.     if ((i % 7) == 3) printf("bn2 div by 7!");
  158.     if (!(i % 11)) printf("bn div by 11!");
  159.     if ((i % 11) == 5) printf("bn2 div by 11!");
  160.     if (!(i % 13)) printf("bn div by 13!");
  161.     if ((i % 13) == 6) printf("bn2 div by 13!");
  162.     i = bnModQ(bn, 7429);    /* 7429 = 17 * 19 * 23 */
  163.     if (!(i % 17)) printf("bn div by 17!");
  164.     if ((i % 17) == 8) printf("bn2 div by 17!");
  165.     if (!(i % 19)) printf("bn div by 19!");
  166.     if ((i % 19) == 9) printf("bn2 div by 19!");
  167.     if (!(i % 23)) printf("bn div by 23!");
  168.     if ((i % 23) == 11) printf("bn2 div by 23!");
  169.     i = bnModQ(bn, 33263);    /* 33263 = 29 * 31 * 37 */
  170.     if (!(i % 29)) printf("bn div by 29!");
  171.     if ((i % 29) == 14) printf("bn2 div by 29!");
  172.     if (!(i % 31)) printf("bn div by 31!");
  173.     if ((i % 31) == 15) printf("bn2 div by 31!");
  174.     if (!(i % 37)) printf("bn div by 37!");
  175.     if ((i % 37) == 18) printf("bn2 div by 37!");
  176. #endif
  177.     /*
  178.      * First, check whether bn is prime.  This uses a fast primality
  179.      * test which usually obviates the need to do one of the
  180.      * confirmation tests later.  See prime.c for a full explanation.
  181.      * We check bn first because it's one bit smaller, saving one
  182.      * modular squaring, and because we might be able to save another
  183.      * when testing it.  (1/4 of the time.)  A small speed hack,
  184.      * but finding big Sophie Germain primes is *slow*.
  185.      */
  186.     if (bnCopy(e, bn) < 0)
  187.         return -1;
  188.     (void)bnSubQ(e, 1);
  189.     l = bnLSWord(e);
  190.  
  191.     j = 1;    /* Where to start in prime array for strong prime tests */
  192.  
  193.     if (l & 7) {
  194.         bnRShift(e, 1);
  195.         if (bnTwoExpMod(a, e, bn) < 0)
  196.             return -1;
  197.         if ((l & 7) == 6) {
  198.             /* bn == 7 mod 8, expect +1 */
  199.             if (bnBits(a) != 1)
  200.                 return 1;    /* Not prime */
  201.             k = 1;
  202.         } else {
  203.             /* bn == 3 or 5 mod 8, expect -1 == bn-1 */
  204.             if (bnAddQ(a, 1) < 0)
  205.                 return -1;
  206.             if (bnCmp(a, bn) != 0)
  207.                 return 1;    /* Not prime */
  208.             k = 1;
  209.             if (l & 4) {
  210.                 /* bn == 5 mod 8, make odd for strong tests */
  211.                 bnRShift(e, 1);
  212.                 k = 2;
  213.             }
  214.         }
  215.     } else {
  216.         /* bn == 1 mod 8, expect 2^((bn-1)/4) == +/-1 mod bn */
  217.         bnRShift(e, 2);
  218.         if (bnTwoExpMod(a, e, bn) < 0)
  219.             return -1;
  220.         if (bnBits(a) == 1) {
  221.             j = 0;    /* Re-do strong prime test to base 2 */
  222.         } else {
  223.             if (bnAddQ(a, 1) < 0)
  224.                 return -1;
  225.             if (bnCmp(a, bn) != 0)
  226.                 return 1;    /* Not prime */
  227.         }
  228.         k = 2 + bnMakeOdd(e);
  229.     }
  230.  
  231.     /*
  232.      * It's prime!  Print a success indicator and check bn2.
  233.      * The success indicator we print is the sign of Jacobi(2,bn2),
  234.      * which is available to us in l.  bn2 = 2*bn + 1.  Since bn is
  235.      * odd, bn2 must be == 3 mod 4, so the options modulo 8 are 3 and 7.
  236.      * 3 if bn == 1 mod 4, 7 if bn == 3 mod 4.  The signs of the
  237.      * Jacobi symbol are - and + in thse two cases.  Set l to be
  238.      * a flag, 1 if the symbol is positive.
  239.      */
  240.     l = (l >> 1) & 1;
  241.     if (f && (err = f(arg, "-+"[l])) < 0)
  242.         return err;
  243.  
  244.     /*
  245.      * Now check bn2.  Since bn2 == 3 mod 4, a strong pseudoprimality
  246.      * test boils down to looking at a^((bn2-1)/2) mod bn and seeing
  247.      * if it's +/-1.  Of course, that exponent is just bn...
  248.      */
  249.     if (bnCopy(bn2, bn) < 0 || bnAdd(bn2, bn) < 0)
  250.         return -1;
  251.     (void)bnAddQ(bn2, 1);    /* Can't overflow */
  252.     if (bnTwoExpMod(a, bn, bn2) < 0)
  253.         return -1;
  254.     if (l) {    /* Expect + */
  255.         if (bnBits(a) != 1)
  256.             return 2;    /* Not prime */
  257.     } else {
  258.         if (bnAddQ(a, 1) < 0)
  259.             return -1;
  260.         if (bnCmp(a, bn2) != 0)
  261.             return 2;    /* Not prime */
  262.     }
  263.  
  264.     /*
  265.      * Success!  We have found a key!  Now go on to confirmation
  266.      * tests...  k is the amount by which e has already been shifted
  267.      * down.  j = 1 unless the test to the base 2 could stand to
  268.      * be re-done, in which case it's 0.
  269.      */
  270.     k += bnMakeOdd(e);
  271.     for (i = j; i < sizeof(primes)/sizeof(*primes); i++) {
  272.         if (f && (err = f(arg, '*')) < 0)
  273.             return err;
  274.  
  275.         /* Check that bn is a strong pseudoprime */
  276.         (void)bnSetQ(a, primes[i]);
  277.         if (bnExpMod(a, a, e, bn) < 0)
  278.             return -1;
  279.  
  280.         if (bnBits(a) != 1) {
  281.             l = k;
  282.             for (;;) {
  283.                 if (bnAddQ(a, 1) < 0)
  284.                     return -1;
  285.                 if (bnCmp(a, bn) == 0)    /* Was result bn-1? */
  286.                     break;    /* Prime */
  287.                 if (!--l)
  288.                     return 2*i+1;    /* Failed, not prime */
  289.                 /* This part is executed once, on average. */
  290.                 (void)bnSubQ(a, 1);    /* Restore a */
  291.                 if (bnSquare(a, a) < 0 || bnMod(a, a, bn) < 0)
  292.                     return -1;
  293.                 if (bnBits(a) == 1)
  294.                     return 2*i+1;    /* Failed, not prime */
  295.             }
  296.         }
  297.  
  298.         /* Okay, that was prime - print success indicator */
  299.         j = bnJacobiQ(primes[i], bn2);
  300.         if (f && (err = f(arg, j < 0 ? '-' : '+')) < 0)
  301.             return err;
  302.         /* If we're re-doing the base 2, we're done. */
  303.         if (!i)
  304.             continue;    /* Already done next part */
  305.  
  306.         /* Check that p^bn == Jacobi(p,bn2) (mod bn2) */
  307.         (void)bnSetQ(a, primes[i]);
  308.         if (bnExpMod(a, a, bn, bn2) < 0)
  309.             return -1;
  310.         /*
  311.          * FIXME:  Actually, we don't need to compute the Jacobi
  312.          * sumbol externally... it never happens that a = +/-1
  313.          * but it's the wrong one.  So we can just look at a and
  314.          * use its sign.  Find a proof somewhere.
  315.          */
  316.         if (j < 0) {
  317.             /* Not a quadratic residue, should have a =  bn2-1 */
  318.             if (bnAddQ(a, 1) < 0)
  319.                 return -1;
  320.             if (bnCmp(a, bn2) != 0)    /* Was result bn2-1? */
  321.                 return 2*i+2;    /* Mismatch -> failed */
  322.         } else {
  323.             /* Quadratic residue, should have a = 1 */
  324.             if (bnBits(a) != 1)
  325.                 return 2*i+2;    /* Mismatch -> failed */
  326.         }
  327.         /* It worked (to the base primes[i]) */
  328.     }
  329.  
  330.     /* Print final success indicator */
  331.     if (f && (err = f(arg, '*')) < 0)
  332.         return err;
  333.     return 0;    /* Prime! */
  334. }
  335.  
  336. /*
  337.  * Modifies the bignum to return the next Sohpie Germain prime >= the
  338.  * input value.  Sohpie Germain primes are number such that p is
  339.  * prime and (p-1)/2 is also prime.
  340.  *
  341.  * Returns >=0 on success or -1 on failure (out of memory).  On success,
  342.  * the return value is the number of modular exponentiations performed
  343.  * (excluding the final confirmations).  This never gives up searching.
  344.  *
  345.  * The FILE *f argument, if non-NULL, has progress indicators written
  346.  * to it.  A dot (.) is written every time a primeality test is failed,
  347.  * a plus (+) or minus (-) when the smaller prime of the pair passes a
  348.  * test, and a star (*) when the larger one does.  Finally, a slash (/)
  349.  * is printed when the sieve was emptied without finding a prime and is
  350.  * being refilled.
  351.  *
  352.  * Apologies to structured programmers for all the GOTOs.
  353.  */
  354. int
  355. germainPrimeGen(struct BigNum *bn2, int (*f)(void *arg, int c), void *arg)
  356. {
  357.     int retval;
  358.     unsigned p, prev;
  359.     struct BigNum a, e, bn;
  360.     int modexps = 0;
  361. #ifdef MSDOS
  362.     unsigned char *sieve;
  363. #else
  364.     unsigned char sieve[SIEVE];
  365. #endif
  366.  
  367. #ifdef MSDOS
  368.     sieve = lbnMemAlloc(SIEVE);
  369.     if (!sieve)
  370.         return -1;
  371. #endif
  372.  
  373.     bnBegin(&a);
  374.     bnBegin(&e);
  375.     bnBegin(&bn);
  376.  
  377.     /*
  378.      * Obviously, the prime we find must be odd.  Further, (p-1)/2
  379.      * must be odd, so p == 3 (mod 3).  Finally, p != 0 (mod 3),
  380.      * and (p-1)/2 != 0 (mod 3), meaning that p != 1 (mod 3).  Thus,
  381.      * p == 11 (mod 12).  Increase bn2 to have this property, and
  382.      * search is steps of 12.  (Actually, we work with the smaller
  383.      * bn, and increase it to 5 mod 6, then search in steps of 6.)
  384.      */
  385.     if (bnCopy(&bn, bn2) < 0)
  386.         goto failed;
  387.     bnRShift(&bn, 1);
  388.     p = bnModQ(&bn, 6);
  389.     if (bnAddQ(&bn, 5-p) < 0)
  390.         goto failed;
  391.  
  392.     for (;;) {
  393.         if (sieveBuild(sieve, SIEVE, &bn, 6, 1) < 0)
  394.             goto failed;
  395.  
  396.         p = prev = 0;
  397.         if (sieve[0] & 1 || (p = sieveSearch(sieve, SIEVE, p)) != 0) {
  398.             do {
  399.                 /*
  400.                  * Adjust bn to have the right value,
  401.                  * incrementing in steps of < 65536.
  402.                  * 10922 = 65536/6.
  403.                  */
  404.                 assert(p >= prev);
  405.                 prev = p-prev;    /* Delta - add 6*prev to bn */
  406. #if SIEVE*8*6 >= 65536
  407.                 while (prev > 10922) {
  408.                     if (bnAddQ(&bn, 6*10922) < 0)
  409.                         goto failed;
  410.                     prev -= 10922;
  411.                 }
  412. #endif
  413.                 if (bnAddQ(&bn, 6*prev) < 0)
  414.                     goto failed;
  415.                 prev = p;
  416.  
  417.                 /* Okay, do the strong tests. */
  418.                 retval = germainPrimeTest(&bn, bn2, &e, &a,
  419.                                           f, arg);
  420.                 if (retval <= 0)
  421.                     goto done;
  422.                 modexps += retval;
  423.                 if (f && (retval = f(arg, '.')) < 0)
  424.                     goto done;
  425.  
  426.                 /* And try again */
  427.                 p = sieveSearch(sieve, SIEVE, p);
  428.             } while (p);
  429.         }
  430.  
  431.         /* Ran out of sieve space - increase bn and keep trying. */
  432. #if SIEVE*8*6 >= 65536
  433.         p = ((SIEVE-1)*8+7) - prev;    /* Number of steps (of 6) */
  434.         while (p >= 10922) {
  435.             if (bnAddQ(&bn, 6*10922) < 0)
  436.                 goto failed;
  437.             p -= 10922;
  438.         }
  439.         if (bnAddQ(&bn, 6*(p+1)) < 0)
  440.             goto failed;
  441. #else
  442.         if (bnAddQ(&bn, SIEVE*8*6 - prev) < 0)
  443.             goto failed;
  444. #endif
  445.         /*
  446.          * Make sure bn2 is up to date for checkpoint/restart
  447.          * operation triggered by calling f with '/'.
  448.          */
  449.         if (bnCopy(bn2, &bn) < 0 || bnAdd(bn2, &bn) < 0)
  450.             goto failed;
  451.         (void)bnAddQ(bn2, 1);    /* Can't overflow */
  452.         if (f && (retval = f(arg, '/')) < 0)
  453.             goto done;
  454.     } /* for (;;) */
  455.  
  456. failed:
  457.     retval = -1;
  458. done:
  459.     bnEnd(&bn);
  460.     bnEnd(&e);
  461.     bnEnd(&a);
  462. #ifdef MSDOS
  463.     lbnMemFree(sieve, SIEVE);
  464. #else
  465.     lbnMemWipe(sieve, sizeof(sieve));
  466. #endif
  467.     return retval < 0 ? retval : modexps;
  468. }
  469.