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