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