![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to keygen.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [PGP] / pgp / src|
Return to keygen.c CVS log | Up to [PGP] / pgp / src |
1.1 root 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:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.