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