home *** CD-ROM | disk | FTP | other *** search
/ PC-Online 1998 February / PCOnline_02_1998.iso / filesbbs / os2 / pgp263.arj / PGP263I.SRC / PGP263II.ZIP / src / genprime.c < prev    next >
C/C++ Source or Header  |  1996-01-02  |  38KB  |  966 lines
  1. /* genprime.c - C source code for generation of large primes
  2.    used by public-key key generation routines.
  3.    First version 17 Mar 87
  4.    Last revised 2 Jun 91 by PRZ
  5.    24 Apr 93 by CP
  6.  
  7.    (c) Copyright 1990-1996 by Philip Zimmermann.  All rights reserved.
  8.    The author assumes no liability for damages resulting from the use
  9.    of this software, even if the damage results from defects in this
  10.    software.  No warranty is expressed or implied.
  11.  
  12.    Note that while most PGP source modules bear Philip Zimmermann's
  13.    copyright notice, many of them have been revised or entirely written
  14.    by contributors who frequently failed to put their names in their
  15.    code.  Code that has been incorporated into PGP from other authors
  16.    was either originally published in the public domain or is used with
  17.    permission from the various authors.
  18.  
  19.    PGP is available for free to the public under certain restrictions.
  20.    See the PGP User's Guide (included in the release package) for
  21.    important information about licensing, patent restrictions on
  22.    certain algorithms, trademarks, copyrights, and export controls.
  23.  
  24.    These functions are for the generation of large prime integers and
  25.    for other functions related to factoring and key generation for 
  26.    many number-theoretic cryptographic algorithms, such as the NIST 
  27.    Digital Signature Standard.
  28.  */
  29.  
  30. #define SHOWPROGRESS
  31.  
  32. /* Define some error status returns for keygen... */
  33. #define NOPRIMEFOUND -14    /* slowtest probably failed */
  34. #define NOSUSPECTS -13        /* fastsieve probably failed */
  35.  
  36.  
  37. #if defined(MSDOS) || defined(WIN32)
  38. #define poll_for_break() {while (kbhit()) getch();}
  39. #endif
  40.  
  41. #ifndef poll_for_break
  42. #define poll_for_break()    /* stub */
  43. #endif
  44.  
  45. #ifdef SHOWPROGRESS
  46. #include <stdio.h>        /* needed for putchar() */
  47. #endif
  48.  
  49. #ifdef MACTC5
  50. extern int  Putchar(int c);
  51. #undef putchar
  52. #define putchar Putchar
  53. #endif
  54.  
  55. #ifdef EMBEDDED            /* compiling for embedded target */
  56. #define _NOMALLOC        /* defined if no malloc is available. */
  57. #endif                /* EMBEDDED */
  58.  
  59. /* Decide whether malloc is available.  Some embedded systems lack it. */
  60. #ifndef _NOMALLOC        /* malloc library routine available */
  61. #include <stdlib.h>        /* ANSI C library - for malloc() and free() */
  62. /* #include <alloc.h> *//* Borland Turbo C has malloc in <alloc.h> */
  63. #endif                /* malloc available */
  64.  
  65. #include "mpilib.h"
  66. #include "genprime.h"
  67. #if (defined(MSDOS) && !defined(__GO32__)) || defined(WIN32)
  68. #include <conio.h>
  69. #endif
  70.  
  71. #include "random.h"
  72.  
  73.  
  74. /* #define STRONGPRIMES *//* if defined, generate "strong" primes for key */
  75. /*
  76.  *"Strong" primes are no longer advantageous, due to the new 
  77.  * elliptical curve method of factoring.  Randomly selected primes 
  78.  * are as good as any.  See "Factoring", by Duncan A. Buell, Journal 
  79.  * of Supercomputing 1 (1987), pages 191-216.
  80.  * This justifies disabling the lengthy search for strong primes.
  81.  *
  82.  * The advice about strong primes in the early RSA literature applies
  83.  * to 256-bit moduli where the attacks were the Pollard rho and P-1
  84.  * factoring algorithms.  Later developments in factoring have entirely
  85.  * supplanted these methods.  The later algorithms are always faster
  86.  * (so we need bigger primes), and don't care about STRONGPRIMES.
  87.  *
  88.  * The early literature was saying that you can get away with small
  89.  * moduli if you choose the primes carefully.  The later developments
  90.  * say you can't get away with small moduli, period.  And it doesn't
  91.  * matter how you choose the primes.
  92.  *
  93.  * It's just taking a heck of a long time for the advice on "strong primes"
  94.  * to disappear from the books.  Authors keep going back to the original
  95.  * documents and repeating what they read there, even though it's out
  96.  * of date.
  97.  */
  98.  
  99. #define BLUM
  100. /* If BLUM is defined, this looks for prines congruent to 3 modulo 4.
  101.    The product of two of these is a Blum integer.  You can uniquely define
  102.    a square root Cmodulo a Blum integer, which leads to some extra
  103.    possibilities for encryption algorithms.  This shrinks the key space by
  104.    2 bits, which is not considered significant.
  105.  */
  106.  
  107. #ifdef STRONGPRIMES
  108.  
  109. static boolean primetest(unitptr p);
  110.     /* Returns TRUE iff p is a prime. */
  111.  
  112. static int mp_sqrt(unitptr quotient, unitptr dividend);
  113.     /* Quotient is returned as the square root of dividend. */
  114.  
  115. #endif
  116.  
  117. static int nextprime(unitptr p);
  118.     /* Find next higher prime starting at p, returning result in p. */
  119.  
  120. static void randombits(unitptr p, short nbits);
  121.     /* Make a random unit array p with nbits of precision. */
  122.  
  123. #ifdef DEBUG
  124. #define DEBUGprintf1(x) printf(x)
  125. #define DEBUGprintf2(x,y) printf(x,y)
  126. #define DEBUGprintf3(x,y,z) printf(x,y,z)
  127. #else
  128. #define DEBUGprintf1(x)
  129. #define DEBUGprintf2(x,y)
  130. #define DEBUGprintf3(x,y,z)
  131. #endif
  132.  
  133.  
  134. /*      primetable is a table of 16-bit prime numbers used for sieving 
  135.    and for other aspects of public-key cryptographic key generation */
  136.  
  137. static word16 primetable[] =
  138. {
  139.     2, 3, 5, 7, 11, 13, 17, 19,
  140.     23, 29, 31, 37, 41, 43, 47, 53,
  141.     59, 61, 67, 71, 73, 79, 83, 89,
  142.     97, 101, 103, 107, 109, 113, 127, 131,
  143.     137, 139, 149, 151, 157, 163, 167, 173,
  144.     179, 181, 191, 193, 197, 199, 211, 223,
  145.     227, 229, 233, 239, 241, 251, 257, 263,
  146.     269, 271, 277, 281, 283, 293, 307, 311,
  147. #ifndef EMBEDDED        /* not embedded, use larger table */
  148.     313, 317, 331, 337, 347, 349, 353, 359,
  149.     367, 373, 379, 383, 389, 397, 401, 409,
  150.     419, 421, 431, 433, 439, 443, 449, 457,
  151.     461, 463, 467, 479, 487, 491, 499, 503,
  152.     509, 521, 523, 541, 547, 557, 563, 569,
  153.     571, 577, 587, 593, 599, 601, 607, 613,
  154.     617, 619, 631, 641, 643, 647, 653, 659,
  155.     661, 673, 677, 683, 691, 701, 709, 719,
  156.     727, 733, 739, 743, 751, 757, 761, 769,
  157.     773, 787, 797, 809, 811, 821, 823, 827,
  158.     829, 839, 853, 857, 859, 863, 877, 881,
  159.     883, 887, 907, 911, 919, 929, 937, 941,
  160.     947, 953, 967, 971, 977, 983, 991, 997,
  161.     1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049,
  162.     1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097,
  163.     1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163,
  164.     1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
  165.     1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283,
  166.     1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321,
  167.     1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423,
  168.     1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
  169.     1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
  170.     1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571,
  171.     1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619,
  172.     1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
  173.     1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
  174.     1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
  175.     1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877,
  176.     1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949,
  177.     1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003,
  178.     2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069,
  179.     2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
  180.     2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203,
  181.     2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
  182.     2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311,
  183.     2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377,
  184.     2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
  185.     2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503,
  186.     2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
  187.     2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657,
  188.     2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693,
  189.     2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
  190.     2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801,
  191.     2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
  192.     2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939,
  193.     2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011,
  194.     3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
  195.     3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167,
  196.     3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221,
  197.     3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301,
  198.     3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347,
  199.     3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
  200.     3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491,
  201.     3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541,
  202.     3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607,
  203.     3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671,
  204.     3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
  205.     3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797,
  206.     3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863,
  207.     3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
  208.     3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003,
  209.     4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
  210.     4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
  211.     4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211,
  212.     4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259,
  213.     4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337,
  214.     4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
  215.     4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481,
  216.     4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547,
  217.     4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621,
  218.     4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673,
  219.     4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
  220.     4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813,
  221.     4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909,
  222.     4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967,
  223.     4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011,
  224.     5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
  225.     5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167,
  226.     5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233,
  227.     5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309,
  228.     5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399,
  229.     5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
  230.     5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507,
  231.     5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573,
  232.     5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653,
  233.     5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711,
  234.     5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
  235.     5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849,
  236.     5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897,
  237.     5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007,
  238.     6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073,
  239.     6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
  240.     6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211,
  241.     6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271,
  242.     6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329,
  243.     6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379,
  244.     6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
  245.     6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563,
  246.     6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637,
  247.     6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701,
  248.     6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779,
  249.     6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
  250.     6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907,
  251.     6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971,
  252.     6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027,
  253.     7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121,
  254.     7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
  255.     7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253,
  256.     7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349,
  257.     7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457,
  258.     7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517,
  259.     7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
  260.     7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621,
  261.     7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691,
  262.     7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757,
  263.     7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853,
  264.     7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
  265.     7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009,
  266.     8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087,
  267.     8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161,
  268.     8167, 8171, 8179, 8191,
  269. #endif                /* not EMBEDDED, use larger table */
  270.     0};            /* null-terminated list, with only one null at end */
  271.  
  272.  
  273.  
  274. #ifdef UNIT8
  275. static word16 bottom16(unitptr r)
  276. /* Called from nextprime and primetest.  Returns low 16 bits of r. */
  277. {
  278.     make_lsbptr(r, (global_precision - ((2 / BYTES_PER_UNIT) - 1)));
  279.     return *(word16 *) r;
  280. }                /* bottom16 */
  281. #else                /* UNIT16 or UNIT32 */
  282. #define bottom16(r) ((word16) lsunit(r))
  283. /* or UNIT32 could mask off lower 16 bits, instead of typecasting it. */
  284. #endif                /* UNIT16 or UNIT32 */
  285.  
  286.  
  287. /*
  288.  * This routine tests p for primality by applying Fermat's theorem:
  289.  * For any x, if ((x**(p-1)) mod p) != 1, then p is not prime.
  290.  * By trying a few values for x, we can determine if p is "probably" prime.
  291.  *
  292.  * Because this test is so slow, it is recommended that p be sieved first
  293.  * to weed out numbers that are obviously not prime.
  294.  *
  295.  * Contrary to what you may have read in the literature, empirical evidence
  296.  * shows this test weeds out a LOT more than 50% of the composite candidates
  297.  * for each trial x.  Each test catches nearly all the composites.
  298.  *
  299.  * Some people have questioned whether four Fermat tests is sufficient.
  300.  * See "Finding Four Million Large Random Primes", by Ronald Rivest,
  301.  * in Advancess in Cryptology: Proceedings of Crypto '91.  He used a
  302.  * small-divisor test similar to PGP's, then a Fermat test to the base 2,
  303.  * and then 8 iterarions of a Miller-Rabin test.  About 718 million random
  304.  * 256-bit integers were generated, 43,741,404 passed the small divisor test,
  305.  * 4,058,000 passed the Fermat test, and all 4,058,000 passed all 8
  306.  * iterations of the Miller-Rabin test, proving their primality beyond most
  307.  * reasonable doubts.  This is strong experimental evidence that the odds
  308.  * of getting a non-prime are less than one in a million (10^-6).
  309.  *
  310.  * He also gives a theoretical argument that the chances of finding a
  311.  * 256-bit non-prime which satisfies one Fermat test to the base 2 is less
  312.  * than 10^-22.  The small divisor test improves this number, and if the
  313.  * numbers are 512 bits (as needed for a 1024-bit key) the odds of failure
  314.  * shrink to about 10^-44.  Thus, he concludes, for practical purposes one
  315.  * Fermat test to the base 2 is sufficient.
  316.  */
  317. static boolean slowtest(unitptr p)
  318. {
  319.     unit x[MAX_UNIT_PRECISION], is_one[MAX_UNIT_PRECISION];
  320.     unit pminus1[MAX_UNIT_PRECISION];
  321.     short i;
  322.  
  323.     mp_move(pminus1, p);
  324.     mp_dec(pminus1);
  325.  
  326.     for (i = 0; i < 4; i++) {    /* Just do a few tests. */
  327.     poll_for_break(); /* polls keyboard, allows ctrl-C to abort program */
  328.     mp_init(x, primetable[i]);    /* Use any old random trial x */
  329.     /* if ((x**(p-1)) mod p) != 1, then p is not prime */
  330.     if (mp_modexp(is_one, x, pminus1, p) < 0)    /* modexp error? */
  331.         return FALSE;    /* error means return not prime status */
  332.     if (testne(is_one, 1))    /* then p is not prime */
  333.         return FALSE;    /* return not prime status */
  334. #ifdef SHOWPROGRESS
  335.     putchar('*');        /* let user see how we are progressing */
  336.     fflush(stdout);
  337. #endif                /* SHOWPROGRESS */
  338.     }
  339.  
  340.     /* If it gets to this point, it's very likely that p is prime */
  341.     mp_burn(x);            /* burn the evidence on the stack... */
  342.     mp_burn(is_one);
  343.     mp_burn(pminus1);
  344.     return TRUE;
  345. }                /* slowtest -- fermattest */
  346.  
  347.  
  348. #ifdef STRONGPRIMES
  349.  
  350. static boolean primetest(unitptr p)
  351. /*
  352.  * Returns TRUE iff p is a prime.
  353.  * If p doesn't pass through the sieve, then p is definitely NOT a prime.
  354.  * If p is small enough for the sieve to prove  primality or not, 
  355.  * and p passes through the sieve, then p is definitely a prime.
  356.  * If p is large and p passes through the sieve and may be a prime,
  357.  * then p is further tested for primality with a slower test.
  358.  */
  359. {
  360.     short i;
  361.     static word16 lastprime = 0;    /* last prime in primetable */
  362.     word16 sqrt_p;    /* to limit sieving past sqrt(p), for small p's */
  363.  
  364.     if (!lastprime) {        /* lastprime still undefined. So define it. */
  365.     /* executes this code only once, then skips it next time */
  366.     for (i = 0; primetable[i]; i++);    /* seek end of primetable */
  367.     lastprime = primetable[i - 1];    /* get last prime in table */
  368.     }
  369.     if (significance(p) <= (2 / BYTES_PER_UNIT))    /* if p <= 16 bits */
  370.     /* p may be in primetable.  Search it. */
  371.     if (bottom16(p) <= lastprime)
  372.         for (i = 0; primetable[i]; i++) {
  373.         /* scan until null-terminator */
  374.         if (primetable[i] == bottom16(p))
  375.             return TRUE;    /* yep, definitely a prime. */
  376.         if (primetable[i] > bottom16(p))    /* we missed. */
  377.             return FALSE;    /* definitely NOT a prime. */
  378.         }        /* We got past the whole primetable without a hit. */
  379.     /* p is bigger than any prime in primetable, so let's sieve. */
  380.     if (!(lsunit(p) & 1))    /* if least significant bit is 0... */
  381.     return FALSE;        /* divisible by 2, not prime */
  382.  
  383.     if (mp_tstminus(p))        /* error if p<0 */
  384.     return FALSE;        /* not prime if p<0 */
  385.  
  386.     /*
  387.      * Optimization for small (32 bits or less) p's.  
  388.      * If p is small, compute sqrt_p = sqrt(p), or else 
  389.      * if p is >32 bits then just set sqrt_p to something 
  390.      * at least as big as the largest primetable entry.
  391.      */
  392.     if (significance(p) <= (4 / BYTES_PER_UNIT)) {    /* if p <= 32 bits */
  393.     unit sqrtp[MAX_UNIT_PRECISION];
  394.     /* Just sieve up to sqrt(p) */
  395.     if (mp_sqrt(sqrtp, p) == 0)    /* 0 means p is a perfect square */
  396.         return FALSE;    /* perfect square is not a prime */
  397.     /* we know that sqrtp <= 16 bits because p <= 32 bits */
  398.     sqrt_p = bottom16(sqrtp);
  399.     } else {
  400.     /* p > 32 bits, so obviate sqrt(p) test. */
  401.     sqrt_p = lastprime;    /* ensures that we do ENTIRE sieve. */
  402.     }
  403.  
  404.     /* p is assumed odd, so begin sieve at 3 */
  405.     for (i = 1; primetable[i]; i++) {
  406.     /* Compute p mod (primetable[i]).  If it divides evenly... */
  407.     if (mp_shortmod(p, primetable[i]) == 0)
  408.         return FALSE;    /* then p is definitely NOT prime */
  409.     if (primetable[i] > sqrt_p)    /* fully sieved p? */
  410.         return TRUE; /* yep, fully passed sieve, definitely a prime. */
  411.     }
  412.     /* It passed the sieve, so p is a suspected prime. */
  413.  
  414.     /*  Now try slow complex primality test on suspected prime. */
  415.     return slowtest(p);        /* returns TRUE or FALSE */
  416. }                /* primetest */
  417.  
  418. #endif
  419.  
  420. /*
  421.  * Used in conjunction with fastsieve.  Builds a table of remainders 
  422.  * relative to the random starting point p, so that fastsieve can 
  423.  * sequentially sieve for suspected primes quickly.  Call buildsieve 
  424.  * once, then call fastsieve for consecutive prime candidates.
  425.  * Note that p must be odd, because the sieve begins at 3. 
  426.  */
  427. static void buildsieve(unitptr p, word16 remainders[])
  428. {
  429.     short i;
  430.     for (i = 1; primetable[i]; i++) {
  431.     remainders[i] = mp_shortmod(p, primetable[i]);
  432.     }
  433. }                /* buildsieve */
  434.  
  435. /*
  436.    Fast prime sieving algorithm by Philip Zimmermann, March 1987.
  437.    This is the fastest algorithm I know of for quickly sieving for 
  438.    large (hundreds of bits in length) "random" probable primes, because 
  439.    it uses only single-precision (16-bit) arithmetic.  Because rigorous 
  440.    prime testing algorithms are very slow, it is recommended that 
  441.    potential prime candidates be quickly passed through this fast 
  442.    sieving algorithm first to weed out numbers that are obviously not 
  443.    prime.
  444.  
  445.    This algorithm is optimized to search sequentially for a large prime 
  446.    from a random starting point.  For generalized nonsequential prime 
  447.    testing, the slower  conventional sieve should be used, as given 
  448.    in primetest(p).
  449.  
  450.    This algorithm requires a fixed table (called primetable) of the 
  451.    first hundred or so small prime numbers. 
  452.    First we select a random odd starting point (p) for our prime 
  453.    search.  Then we build a table of 16-bit remainders calculated 
  454.    from that initial p.  This table of 16-bit remainders is exactly 
  455.    the same length as the table of small 16-bit primes.  Each 
  456.    remainders table entry contains the remainder of p divided by the 
  457.    corresponding primetable entry.  Then we begin sequentially testing 
  458.    all odd integers, starting from the initial odd random p.  The 
  459.    distance we have searched from the huge random starting point p is 
  460.    a small 16-bit number, pdelta.  If pdelta plus a remainders table 
  461.    entry is evenly divisible by the corresponding primetable entry, 
  462.    then p+pdelta is factorable by that primetable entry, which means 
  463.    p+pdelta is not prime.
  464.  */
  465.  
  466. /*      Fastsieve is used for searching sequentially from a random starting
  467.    point for a suspected prime.  Requires that buildsieve be called 
  468.    first, to build a table of remainders relative to the random starting 
  469.    point p.  
  470.    Returns true iff pdelta passes through the sieve and thus p+pdelta 
  471.    may be a prime.  Note that p must be odd, because the sieve begins 
  472.    at 3.
  473.  */
  474. static boolean fastsieve(word16 pdelta, word16 remainders[])
  475. {
  476.     short i;
  477.     for (i = 1; primetable[i]; i++) {
  478.     /*
  479.      * If pdelta plus a remainders table entry is evenly 
  480.      * divisible by the corresponding primetable entry,
  481.      * then p+pdelta is factorable by that primetable entry, 
  482.      * which means p+pdelta is not prime.
  483.      */
  484.     if ((pdelta + remainders[i]) % primetable[i] == 0)
  485.         return FALSE;    /* then p+pdelta is not prime */
  486.     }
  487.     /* It passed the sieve.  It is now a suspected prime. */
  488.     return TRUE;
  489. }                /* fastsieve */
  490.  
  491.  
  492. #define numberof(x) (sizeof(x)/sizeof(x[0]))    /* number of table entries */
  493.  
  494.  
  495. static int nextprime(unitptr p)
  496. /*
  497.  * Find next higher prime starting at p, returning result in p. 
  498.  * Uses fast prime sieving algorithm to search sequentially.
  499.  * Returns 0 for normal completion status, < 0 for failure status.
  500.  */
  501. {
  502.     word16 pdelta, range;
  503.     short oldprecision;
  504.     short i, suspects;
  505.  
  506.     /* start search at candidate p */
  507.     mp_inc(p);    /* remember, it's the NEXT prime from p, noninclusive. */
  508.     if (significance(p) <= 1) {
  509.     /*
  510.      * p might be smaller than the largest prime in primetable.
  511.      * We can't sieve for primes that are already in primetable.
  512.      * We will have to directly search the table.
  513.      */
  514.     /* scan until null-terminator */
  515.     for (i = 0; primetable[i]; i++) {
  516.         if (primetable[i] >= lsunit(p)) {
  517.         mp_init(p, primetable[i]);
  518.         return 0;    /* return next higher prime from primetable */
  519.         }
  520.     }        /* We got past the whole primetable without a hit. */
  521.     }          /* p is bigger than any prime in primetable, so let's sieve. */
  522.     if (mp_tstminus(p)) {    /* error if p<0 */
  523.     mp_init(p, 2);        /* next prime >0 is 2 */
  524.     return 0;        /* normal completion status */
  525.     }
  526. #ifndef BLUM
  527.     lsunit(p) |= 1;        /* set candidate's lsb - make it odd */
  528. #else
  529.     lsunit(p) |= 3;        /* Make candidate ==3 mod 4 */
  530. #endif
  531.  
  532.     /* Adjust the global_precision downward to the optimum size for p... */
  533.     oldprecision = global_precision;    /* save global_precision */
  534.     /* We will need 2-3 extra bits of precision for the falsekeytest. */
  535.     set_precision(bits2units(countbits(p) + 4 + SLOP_BITS));
  536.     /* Rescale p to global_precision we just defined */
  537.     rescale(p, oldprecision, global_precision);
  538.  
  539.     {
  540. #ifdef _NOMALLOC /* No malloc and free functions available.  Use stack. */
  541.     word16 remainders[numberof(primetable)];
  542. #else            /* malloc available, we can conserve stack space. */
  543.     word16 *remainders;
  544.     /* Allocate some memory for the table of remainders: */
  545.     remainders = (word16 *) malloc(sizeof(primetable));
  546. #endif                /* malloc available */
  547.  
  548.     /* Build remainders table relative to initial p: */
  549.     buildsieve(p, remainders);
  550.     pdelta = 0;        /* offset from initial random p */
  551.     /* Sieve preparation complete.  Now for some fast fast sieving... */
  552.     /* slowtest will not be called unless fastsieve is true */
  553.  
  554.     /* range is how far to search before giving up. */
  555. #ifndef BLUM
  556.     range = 4 * units2bits(global_precision);
  557. #else
  558.     /* Twice as many because step size is twice as large, */
  559.     range = 8 * units2bits(global_precision);
  560. #endif
  561.     suspects = 0;    /* number of suspected primes and slowtest trials */
  562.     for (;;) {
  563.         /* equivalent to:  if (primetest(p)) break; */
  564.         if (fastsieve(pdelta, remainders)) { /* found suspected prime */
  565.         suspects++;    /* tally for statistical purposes */
  566. #ifdef SHOWPROGRESS
  567.         putchar('.');    /* let user see how we are progressing */
  568.         fflush(stdout);
  569. #endif                /* SHOWPROGRESS */
  570.         if (slowtest(p))
  571.             break;    /* found a prime */
  572.         }
  573. #ifndef BLUM
  574.         pdelta += 2;    /* try next odd number */
  575. #else
  576.         pdelta += 4;
  577.         mp_inc(p);
  578.         mp_inc(p);
  579. #endif
  580.         mp_inc(p);
  581.         mp_inc(p);
  582.  
  583.         if (pdelta > range)    /* searched too many candidates? */
  584.         break;    /* something must be wrong--bail out of search */
  585.  
  586.     }            /* while (TRUE) */
  587.  
  588. #ifdef SHOWPROGRESS
  589.     putchar(' ');        /* let user see how we are progressing */
  590. #endif                /* SHOWPROGRESS */
  591.  
  592.     for (i = 0; primetable[i]; i++)    /* scan until null-terminator */
  593.         remainders[i] = 0;    /* don't leave remainders exposed in RAM */
  594. #ifndef _NOMALLOC
  595.     free(remainders);    /* free allocated memory */
  596. #endif                /* not _NOMALLOC */
  597.     }
  598.  
  599.     set_precision(oldprecision);    /* restore precision */
  600.  
  601.     if (pdelta > range) {    /* searched too many candidates? */
  602.     if (suspects < 1)    /* unreasonable to have found no suspects */
  603.         return NOSUSPECTS;    /* fastsieve failed, probably */
  604.     return NOPRIMEFOUND;    /* return error status */
  605.     }
  606.     return 0;            /* return normal completion status */
  607.  
  608. }                /* nextprime */
  609.  
  610.  
  611. /* We will need a series of truly random bits for key generation.
  612.    In most implementations, our random number supply is derived from
  613.    random keyboard delays rather than a hardware random number
  614.    chip.  So we will have to ensure we have a large enough pool of
  615.    accumulated random numbers from the keyboard.  trueRandAccum()
  616.    performs this operation.  
  617.  */
  618.  
  619. /* Fills 1 unit with random bytes, and returns unit. */
  620. static unit randomunit(void)
  621. {
  622.     unit u = 0;
  623.     byte i;
  624.     i = BYTES_PER_UNIT;
  625.     do
  626.     u = (u << 8) + trueRandByte();
  627.     while (--i != 0);
  628.     return u;
  629. }                /* randomunit */
  630.  
  631. /*
  632.  * Make a random unit array p with nbits of precision.  Used mainly to 
  633.  * generate large random numbers to search for primes.
  634.  */
  635. static void randombits(unitptr p, short nbits)
  636. {
  637.     mp_init(p, 0);
  638.     make_lsbptr(p, global_precision);
  639.  
  640.     /* Add whole units of randomness */
  641.     while (nbits >= UNITSIZE) {
  642.     *post_higherunit(p) = randomunit();
  643.     nbits -= UNITSIZE;
  644.     }
  645.  
  646.     /* Add most-significant partial unit (if any) */
  647.     if (nbits)
  648.     *p = randomunit() & (power_of_2(nbits) - 1);
  649. }                /* randombits */
  650.  
  651. /*
  652.  * Makes a "random" prime p with nbits significant bits of precision.
  653.  * Since these primes are used to compute a modulus of a guaranteed 
  654.  * length, the top 2 bits of the prime are set to 1, so that the
  655.  * product of 2 primes (the modulus) is of a deterministic length.
  656.  * Returns 0 for normal completion status, < 0 for failure status.
  657.  */
  658. int randomprime(unitptr p, short nbits)
  659. {
  660.     DEBUGprintf2("\nGenerating a %d-bit random prime. ", nbits);
  661.     /* Get an initial random candidate p to start search. */
  662.     randombits(p, nbits - 2);    /* 2 less random bits for nonrandom top bits */
  663.     /* To guarantee exactly nbits of significance, set the top 2 bits to 1 */
  664.     mp_setbit(p, nbits - 1);    /* highest bit is nonrandom */
  665.     mp_setbit(p, nbits - 2);    /* next highest bit is also nonrandom */
  666.     return nextprime(p);    /* search for next higher prime
  667.                    from starting point p */
  668. }                /* randomprime */
  669.  
  670.  
  671. #ifdef STRONGPRIMES        /* generate "strong" primes for keys */
  672.  
  673. #define log_1stprime 6        /* log base 2 of firstprime */
  674.  
  675. /* 1st primetable entry used by tryprime */
  676. #define firstprime (1<<log_1stprime)
  677.  
  678. /* This routine attempts to generate a prime p such that p-1 has prime p1
  679.    as its largest factor.  Prime p will have no more than maxbits bits of
  680.    significance.  Prime p1 must be less than maxbits-log_1stprime in length.  
  681.    This routine is called only from goodprime.
  682.  */
  683. static boolean tryprime(unitptr p, unitptr p1, short maxbits)
  684. {
  685.     int i;
  686.     unit i2[MAX_UNIT_PRECISION];
  687.     /* Generate p such that p = (i*2*p1)+1, for i=1,2,3,5,7,11,13,17...
  688.        and test p for primality for each small prime i.
  689.        It's better to start i at firstprime rather than at 1,
  690.        because then p grows slower in significance.
  691.        Start looking for small primes that are > firstprime...
  692.      */
  693.     if ((countbits(p1) + log_1stprime) >= maxbits) {
  694.     DEBUGprintf1("\007[Error: overconstrained prime]");
  695.     return FALSE;        /* failed to make a good prime */
  696.     }
  697.     for (i = 0; primetable[i]; i++) {
  698.     if (primetable[i] < firstprime)
  699.         continue;
  700.     /* note that mp_init doesn't extend sign bit for >32767 */
  701.     mp_init(i2, primetable[i] << 1);
  702.     mp_mult(p, p1, i2);
  703.     mp_inc(p);
  704.     if (countbits(p) > maxbits)
  705.         break;
  706.     DEBUGprintf1(".");
  707.     if (primetest(p))
  708.         return TRUE;
  709.     }
  710.     return FALSE;        /* failed to make a good prime */
  711. }                /* tryprime */
  712.  
  713.  
  714. /*
  715.  * Make a "strong" prime p with at most maxbits and at least minbits of 
  716.  * significant bits of precision.  This algorithm is called to generate
  717.  * a high-quality prime p for key generation purposes.  It must have 
  718.  * special characteristics for making a modulus n that is hard to factor.
  719.  * Returns 0 for normal completion status, < 0 for failure status.
  720.  */
  721. int goodprime(unitptr p, short maxbits, short minbits)
  722. {
  723.     unit p1[MAX_UNIT_PRECISION];
  724.     short oldprecision, midbits;
  725.     int status;
  726.  
  727.     mp_init(p, 0);
  728.     /* Adjust the global_precision downward to the optimum size for p... */
  729.     oldprecision = global_precision;    /* save global_precision */
  730.     /* We will need 2-3 extra bits of precision for the falsekeytest. */
  731.     set_precision(bits2units(maxbits + 4 + SLOP_BITS));
  732.     /* rescale p to global_precision we just defined */
  733.     rescale(p, oldprecision, global_precision);
  734.  
  735.     minbits -= 2 * log_1stprime;    /* length of p" */
  736.     midbits = (maxbits + minbits) / 2;    /* length of p' */
  737.     DEBUGprintf3("\nGenerating a %d-%d bit refined prime. ",
  738.          minbits + 2 * log_1stprime, maxbits);
  739.     do {
  740.     do {
  741.         status = randomprime(p, minbits - 1);
  742.         if (status < 0)
  743.         return status;    /* failed to find a random prime */
  744.         DEBUGprintf2("\n(p\042=%d bits)", countbits(p));
  745.     } while (!tryprime(p1, p, midbits));
  746.     DEBUGprintf2("(p'=%d bits)", countbits(p1));
  747.     } while (!tryprime(p, p1, maxbits));
  748.     DEBUGprintf2("\n\007(p=%d bits) ", countbits(p));
  749.     mp_burn(p1);        /* burn the evidence on the stack */
  750.     set_precision(oldprecision);    /* restore precision */
  751.     return 0;            /* normal completion status */
  752. }                /* goodprime */
  753.  
  754. #endif                /* STRONGPRIMES */
  755.  
  756.  
  757. #define iplus1  ( i==2 ? 0 : i+1 )    /* used by Euclid algorithms */
  758. #define iminus1 ( i==0 ? 2 : i-1 )    /* used by Euclid algorithms */
  759.  
  760. /* Computes greatest common divisor via Euclid's algorithm. */
  761. void mp_gcd(unitptr result, unitptr a, unitptr n)
  762. {
  763.     short i;
  764.     unit gcopies[3][MAX_UNIT_PRECISION];
  765. #define g(i) (  &(gcopies[i][0])  )
  766.     mp_move(g(0), n);
  767.     mp_move(g(1), a);
  768.  
  769.     i = 1;
  770.     while (testne(g(i), 0)) {
  771.     mp_mod(g(iplus1), g(iminus1), g(i));
  772.     i = iplus1;
  773.     }
  774.     mp_move(result, g(iminus1));
  775.     mp_burn(g(iminus1));    /* burn the evidence on the stack... */
  776.     mp_burn(g(iplus1));
  777. #undef g
  778. }                /* mp_gcd */
  779.  
  780. /*
  781.  * Euclid's algorithm extended to compute multiplicative inverse.
  782.  * Computes x such that a*x mod n = 1, where 0<a<n
  783.  *
  784.  * The variable u is unnecessary for the algorithm, but is 
  785.  * included in comments for mathematical clarity. 
  786.  */
  787. void mp_inv(unitptr x, unitptr a, unitptr n)
  788. {
  789.     short i;
  790.     unit y[MAX_UNIT_PRECISION], temp[MAX_UNIT_PRECISION];
  791.     unit gcopies[3][MAX_UNIT_PRECISION], vcopies[3][MAX_UNIT_PRECISION];
  792. #define g(i) (  &(gcopies[i][0])  )
  793. #define v(i) (  &(vcopies[i][0])  )
  794. /*      unit ucopies[3][MAX_UNIT_PRECISION]; */
  795. /* #define u(i) (  &(ucopies[i][0])  ) */
  796.     mp_move(g(0), n);
  797.     mp_move(g(1), a);
  798. /*      mp_init(u(0),1); mp_init(u(1),0); */
  799.     mp_init(v(0), 0);
  800.     mp_init(v(1), 1);
  801.     i = 1;
  802.     while (testne(g(i), 0)) {
  803.     /* we know that at this point,  g(i) = u(i)*n + v(i)*a  */
  804.     mp_udiv(g(iplus1), y, g(iminus1), g(i));
  805.     mp_mult(temp, y, v(i));
  806.     mp_move(v(iplus1), v(iminus1));
  807.     mp_sub(v(iplus1), temp);
  808. /*      mp_mult(temp,y,u(i)); mp_move(u(iplus1),u(iminus1));
  809.     mp_sub(u(iplus1),temp); */
  810.     i = iplus1;
  811.     }
  812.     mp_move(x, v(iminus1));
  813.     if (mp_tstminus(x))
  814.     mp_add(x, n);
  815.     mp_burn(g(iminus1));    /* burn the evidence on the stack... */
  816.     mp_burn(g(iplus1));
  817.     mp_burn(v(0));
  818.     mp_burn(v(1));
  819.     mp_burn(v(2));
  820.     mp_burn(y);
  821.     mp_burn(temp);
  822. #undef g
  823. #undef v
  824. }                /* mp_inv */
  825.  
  826. #ifdef STRONGPRIMES
  827.  
  828. /*      mp_sqrt - returns square root of a number.
  829.    returns -1 for error, 0 for perfect square, 1 for not perfect square.
  830.    Not used by any RSA-related functions.       Used by factoring algorithms.
  831.    This version needs optimization.
  832.    by Charles W. Merritt  July 15, 1989, refined by PRZ.
  833.  
  834.    These are notes on computing the square root the manual old-fashioned 
  835.    way.  This is the basis of the fast sqrt algorithm mp_sqrt below:
  836.  
  837.    1)   Separate the number into groups (periods) of two digits each,
  838.    beginning with units or at the decimal point.
  839.    2)   Find the greatest perfect square in the left hand period & write 
  840.    its  square root as the first figure of the required root.  Subtract
  841.    the square of this number from the left hand period.  Annex to the
  842.    remainder the next group so as to form a dividend.
  843.    3)   Double the root already found and write it as a partial divisor at 
  844.    the left of the new dividend.  Annex one zero digit, making a trial 
  845.    divisor, and divide the new dividend by the trial divisor.
  846.    4)   Write the quotient in the root as the trial term and also substitute 
  847.    this quotient for the annexed zero digit in the partial divisor, 
  848.    making the latter complete.
  849.    5)   Multiply the complete divisor by the figure just obtained and, 
  850.    if possible, subtract the product from the last remainder.
  851.    If this product is too large, the trial term of the quotient
  852.    must be replaced by the next smaller number and the operations
  853.    preformed as before.
  854.    (IN BINARY, OUR TRIAL TERM IS ALWAYS 1 AND WE USE IT OR DO NOTHING.)
  855.    6)   Proceed in this manner until all periods are used.
  856.    If there is still a remainder, it's not a perfect square.
  857.  */
  858.  
  859. /* Quotient is returned as the square root of dividend. */
  860. static int mp_sqrt(unitptr quotient, unitptr dividend)
  861. {
  862.     register short next2bits;    /* "period", or group of 2 bits of dividend */
  863.     register unit dvdbitmask, qbitmask;
  864.     unit remainder[MAX_UNIT_PRECISION];
  865.     unit rjq[MAX_UNIT_PRECISION], divisor[MAX_UNIT_PRECISION];
  866.     unsigned int qbits, qprec, dvdbits, dprec, oldprecision;
  867.     int notperfect;
  868.  
  869.     mp_init(quotient, 0);
  870.     if (mp_tstminus(dividend)) {    /* if dividend<0, return error */
  871.     mp_dec(quotient);    /* quotient = -1 */
  872.     return -1;
  873.     }
  874.     /* normalize and compute number of bits in dividend first */
  875.     init_bitsniffer(dividend, dvdbitmask, dprec, dvdbits);
  876.     /* init_bitsniffer returns a 0 if dvdbits is 0 */
  877.     if (dvdbits == 1) {
  878.     mp_init(quotient, 1);    /* square root of 1 is 1 */
  879.     return 0;
  880.     }
  881.     /* rescale quotient to half the precision of dividend */
  882.     qbits = (dvdbits + 1) >> 1;
  883.     qprec = bits2units(qbits);
  884.     rescale(quotient, global_precision, qprec);
  885.     make_msbptr(quotient, qprec);
  886.     qbitmask = power_of_2((qbits - 1) & (UNITSIZE - 1));
  887.  
  888.     /*
  889.      * Set smallest optimum precision for this square root.
  890.      * The low-level primitives are affected by the call to set_precision.
  891.      * Even though the dividend precision is bigger than the precision
  892.      * we will use, no low-level primitives will be used on the dividend.
  893.      * They will be used on the quotient, remainder, and rjq, which are
  894.      * smaller precision.
  895.      */
  896.     oldprecision = global_precision;    /* save global_precision */
  897.     set_precision(bits2units(qbits + 3));    /* 3 bits of precision slop */
  898.  
  899.     /* special case: sqrt of 1st 2 (binary) digits of dividend
  900.        is 1st (binary) digit of quotient.  This is always 1. */
  901.     stuff_bit(quotient, qbitmask);
  902.     bump_bitsniffer(quotient, qbitmask);
  903.     mp_init(rjq, 1);        /* rjq is Right Justified Quotient */
  904.  
  905.     if (!(dvdbits & 1)) {
  906.     /* even number of bits in dividend */
  907.     next2bits = 2;
  908.     bump_bitsniffer(dividend, dvdbitmask);
  909.     dvdbits--;
  910.     if (sniff_bit(dividend, dvdbitmask))
  911.         next2bits++;
  912.     bump_bitsniffer(dividend, dvdbitmask);
  913.     dvdbits--;
  914.     } else {
  915.     /* odd number of bits in dividend */
  916.     next2bits = 1;
  917.     bump_bitsniffer(dividend, dvdbitmask);
  918.     dvdbits--;
  919.     }
  920.  
  921.     mp_init(remainder, next2bits - 1);
  922.  
  923.     /* dvdbits is guaranteed to be even at this point */
  924.  
  925.     while (dvdbits) {
  926.     next2bits = 0;
  927.     if (sniff_bit(dividend, dvdbitmask))
  928.         next2bits = 2;
  929.     bump_bitsniffer(dividend, dvdbitmask);
  930.     dvdbits--;
  931.     if (sniff_bit(dividend, dvdbitmask))
  932.         next2bits++;
  933.     bump_bitsniffer(dividend, dvdbitmask);
  934.     dvdbits--;
  935.     mp_rotate_left(remainder, (boolean) ((next2bits & 2) != 0));
  936.     mp_rotate_left(remainder, (boolean) ((next2bits & 1) != 0));
  937.  
  938.     /*
  939.      * "divisor" is trial divisor, complete divisor is 4*rjq 
  940.      * or 4*rjq+1.
  941.      * Subtract divisor times its last digit from remainder.
  942.      * If divisor ends in 1, remainder -= divisor*1,
  943.      * or if divisor ends in 0, remainder -= divisor*0 (do nothing).
  944.      * Last digit of divisor inflates divisor as large as possible
  945.      * yet still subtractable from remainder.
  946.      */
  947.     mp_move(divisor, rjq);    /* divisor = 4*rjq+1 */
  948.     mp_rotate_left(divisor, 0);
  949.     mp_rotate_left(divisor, 1);
  950.     if (mp_compare(remainder, divisor) >= 0) {
  951.         mp_sub(remainder, divisor);
  952.         stuff_bit(quotient, qbitmask);
  953.         mp_rotate_left(rjq, 1);
  954.     } else {
  955.         mp_rotate_left(rjq, 0);
  956.     }
  957.     bump_bitsniffer(quotient, qbitmask);
  958.     }
  959.     notperfect = testne(remainder, 0);    /* not a perfect square? */
  960.     set_precision(oldprecision);    /* restore original precision */
  961.     return notperfect;        /* normal return */
  962. }                /* mp_sqrt */
  963. #endif
  964.  
  965. /*------------------- End of genprime.c -----------------------------*/
  966.