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