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