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