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