1 /*
2 * mpprime.c
3 *
4 * Utilities for finding and working with prime and pseudo-prime
5 * integers
6 *
7 * This Source Code Form is subject to the terms of the Mozilla Public
8 * License, v. 2.0. If a copy of the MPL was not distributed with this
9 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
10
11 #include "mpi-priv.h"
12 #include "mpprime.h"
13 #include "mplogic.h"
14 #include <stdlib.h>
15 #include <string.h>
16
17 #define SMALL_TABLE 0 /* determines size of hard-wired prime table */
18
19 #define RANDOM() rand()
20
21 #include "primes.c" /* pull in the prime digit table */
22
23 /*
24 Test if any of a given vector of digits divides a. If not, MP_NO
25 is returned; otherwise, MP_YES is returned and 'which' is set to
26 the index of the integer in the vector which divided a.
27 */
28 mp_err s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which);
29
30 /* {{{ mpp_divis(a, b) */
31
32 /*
33 mpp_divis(a, b)
34
35 Returns MP_YES if a is divisible by b, or MP_NO if it is not.
36 */
37
38 mp_err
mpp_divis(mp_int * a,mp_int * b)39 mpp_divis(mp_int *a, mp_int *b)
40 {
41 mp_err res;
42 mp_int rem;
43
44 if ((res = mp_init(&rem)) != MP_OKAY)
45 return res;
46
47 if ((res = mp_mod(a, b, &rem)) != MP_OKAY)
48 goto CLEANUP;
49
50 if (mp_cmp_z(&rem) == 0)
51 res = MP_YES;
52 else
53 res = MP_NO;
54
55 CLEANUP:
56 mp_clear(&rem);
57 return res;
58
59 } /* end mpp_divis() */
60
61 /* }}} */
62
63 /* {{{ mpp_divis_d(a, d) */
64
65 /*
66 mpp_divis_d(a, d)
67
68 Return MP_YES if a is divisible by d, or MP_NO if it is not.
69 */
70
71 mp_err
mpp_divis_d(mp_int * a,mp_digit d)72 mpp_divis_d(mp_int *a, mp_digit d)
73 {
74 mp_err res;
75 mp_digit rem;
76
77 ARGCHK(a != NULL, MP_BADARG);
78
79 if (d == 0)
80 return MP_NO;
81
82 if ((res = mp_mod_d(a, d, &rem)) != MP_OKAY)
83 return res;
84
85 if (rem == 0)
86 return MP_YES;
87 else
88 return MP_NO;
89
90 } /* end mpp_divis_d() */
91
92 /* }}} */
93
94 /* {{{ mpp_random(a) */
95
96 /*
97 mpp_random(a)
98
99 Assigns a random value to a. This value is generated using the
100 standard C library's rand() function, so it should not be used for
101 cryptographic purposes, but it should be fine for primality testing,
102 since all we really care about there is good statistical properties.
103
104 As many digits as a currently has are filled with random digits.
105 */
106
107 mp_err
mpp_random(mp_int * a)108 mpp_random(mp_int *a)
109
110 {
111 mp_digit next = 0;
112 unsigned int ix, jx;
113
114 ARGCHK(a != NULL, MP_BADARG);
115
116 for (ix = 0; ix < USED(a); ix++) {
117 for (jx = 0; jx < sizeof(mp_digit); jx++) {
118 next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX);
119 }
120 DIGIT(a, ix) = next;
121 }
122
123 return MP_OKAY;
124
125 } /* end mpp_random() */
126
127 /* }}} */
128
129 /* {{{ mpp_random_size(a, prec) */
130
131 mp_err
mpp_random_size(mp_int * a,mp_size prec)132 mpp_random_size(mp_int *a, mp_size prec)
133 {
134 mp_err res;
135
136 ARGCHK(a != NULL && prec > 0, MP_BADARG);
137
138 if ((res = s_mp_pad(a, prec)) != MP_OKAY)
139 return res;
140
141 return mpp_random(a);
142
143 } /* end mpp_random_size() */
144
145 /* }}} */
146
147 /* {{{ mpp_divis_vector(a, vec, size, which) */
148
149 /*
150 mpp_divis_vector(a, vec, size, which)
151
152 Determines if a is divisible by any of the 'size' digits in vec.
153 Returns MP_YES and sets 'which' to the index of the offending digit,
154 if it is; returns MP_NO if it is not.
155 */
156
157 mp_err
mpp_divis_vector(mp_int * a,const mp_digit * vec,int size,int * which)158 mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which)
159 {
160 ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG);
161
162 return s_mpp_divp(a, vec, size, which);
163
164 } /* end mpp_divis_vector() */
165
166 /* }}} */
167
168 /* {{{ mpp_divis_primes(a, np) */
169
170 /*
171 mpp_divis_primes(a, np)
172
173 Test whether a is divisible by any of the first 'np' primes. If it
174 is, returns MP_YES and sets *np to the value of the digit that did
175 it. If not, returns MP_NO.
176 */
177 mp_err
mpp_divis_primes(mp_int * a,mp_digit * np)178 mpp_divis_primes(mp_int *a, mp_digit *np)
179 {
180 int size, which;
181 mp_err res;
182
183 ARGCHK(a != NULL && np != NULL, MP_BADARG);
184
185 size = (int)*np;
186 if (size > prime_tab_size)
187 size = prime_tab_size;
188
189 res = mpp_divis_vector(a, prime_tab, size, &which);
190 if (res == MP_YES)
191 *np = prime_tab[which];
192
193 return res;
194
195 } /* end mpp_divis_primes() */
196
197 /* }}} */
198
199 /* {{{ mpp_fermat(a, w) */
200
201 /*
202 Using w as a witness, try pseudo-primality testing based on Fermat's
203 little theorem. If a is prime, and (w, a) = 1, then w^a == w (mod
204 a). So, we compute z = w^a (mod a) and compare z to w; if they are
205 equal, the test passes and we return MP_YES. Otherwise, we return
206 MP_NO.
207 */
208 mp_err
mpp_fermat(mp_int * a,mp_digit w)209 mpp_fermat(mp_int *a, mp_digit w)
210 {
211 mp_int base, test;
212 mp_err res;
213
214 if ((res = mp_init(&base)) != MP_OKAY)
215 return res;
216
217 mp_set(&base, w);
218
219 if ((res = mp_init(&test)) != MP_OKAY)
220 goto TEST;
221
222 /* Compute test = base^a (mod a) */
223 if ((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY)
224 goto CLEANUP;
225
226 if (mp_cmp(&base, &test) == 0)
227 res = MP_YES;
228 else
229 res = MP_NO;
230
231 CLEANUP:
232 mp_clear(&test);
233 TEST:
234 mp_clear(&base);
235
236 return res;
237
238 } /* end mpp_fermat() */
239
240 /* }}} */
241
242 /*
243 Perform the fermat test on each of the primes in a list until
244 a) one of them shows a is not prime, or
245 b) the list is exhausted.
246 Returns: MP_YES if it passes tests.
247 MP_NO if fermat test reveals it is composite
248 Some MP error code if some other error occurs.
249 */
250 mp_err
mpp_fermat_list(mp_int * a,const mp_digit * primes,mp_size nPrimes)251 mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes)
252 {
253 mp_err rv = MP_YES;
254
255 while (nPrimes-- > 0 && rv == MP_YES) {
256 rv = mpp_fermat(a, *primes++);
257 }
258 return rv;
259 }
260
261 /* {{{ mpp_pprime(a, nt) */
262
263 /*
264 mpp_pprime(a, nt)
265
266 Performs nt iteration of the Miller-Rabin probabilistic primality
267 test on a. Returns MP_YES if the tests pass, MP_NO if one fails.
268 If MP_NO is returned, the number is definitely composite. If MP_YES
269 is returned, it is probably prime (but that is not guaranteed).
270 */
271
272 mp_err
mpp_pprime(mp_int * a,int nt)273 mpp_pprime(mp_int *a, int nt)
274 {
275 mp_err res;
276 mp_int x, amo, m, z; /* "amo" = "a minus one" */
277 int iter;
278 unsigned int jx;
279 mp_size b;
280
281 ARGCHK(a != NULL, MP_BADARG);
282
283 MP_DIGITS(&x) = 0;
284 MP_DIGITS(&amo) = 0;
285 MP_DIGITS(&m) = 0;
286 MP_DIGITS(&z) = 0;
287
288 /* Initialize temporaries... */
289 MP_CHECKOK(mp_init(&amo));
290 /* Compute amo = a - 1 for what follows... */
291 MP_CHECKOK(mp_sub_d(a, 1, &amo));
292
293 b = mp_trailing_zeros(&amo);
294 if (!b) { /* a was even ? */
295 res = MP_NO;
296 goto CLEANUP;
297 }
298
299 MP_CHECKOK(mp_init_size(&x, MP_USED(a)));
300 MP_CHECKOK(mp_init(&z));
301 MP_CHECKOK(mp_init(&m));
302 MP_CHECKOK(mp_div_2d(&amo, b, &m, 0));
303
304 /* Do the test nt times... */
305 for (iter = 0; iter < nt; iter++) {
306
307 /* Choose a random value for 1 < x < a */
308 MP_CHECKOK(s_mp_pad(&x, USED(a)));
309 mpp_random(&x);
310 MP_CHECKOK(mp_mod(&x, a, &x));
311 if (mp_cmp_d(&x, 1) <= 0) {
312 iter--; /* don't count this iteration */
313 continue; /* choose a new x */
314 }
315
316 /* Compute z = (x ** m) mod a */
317 MP_CHECKOK(mp_exptmod(&x, &m, a, &z));
318
319 if (mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) {
320 res = MP_YES;
321 continue;
322 }
323
324 res = MP_NO; /* just in case the following for loop never executes. */
325 for (jx = 1; jx < b; jx++) {
326 /* z = z^2 (mod a) */
327 MP_CHECKOK(mp_sqrmod(&z, a, &z));
328 res = MP_NO; /* previous line set res to MP_YES */
329
330 if (mp_cmp_d(&z, 1) == 0) {
331 break;
332 }
333 if (mp_cmp(&z, &amo) == 0) {
334 res = MP_YES;
335 break;
336 }
337 } /* end testing loop */
338
339 /* If the test passes, we will continue iterating, but a failed
340 test means the candidate is definitely NOT prime, so we will
341 immediately break out of this loop
342 */
343 if (res == MP_NO)
344 break;
345
346 } /* end iterations loop */
347
348 CLEANUP:
349 mp_clear(&m);
350 mp_clear(&z);
351 mp_clear(&x);
352 mp_clear(&amo);
353 return res;
354
355 } /* end mpp_pprime() */
356
357 /* }}} */
358
359 /* Produce table of composites from list of primes and trial value.
360 ** trial must be odd. List of primes must not include 2.
361 ** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest
362 ** prime in list of primes. After this function is finished,
363 ** if sieve[i] is non-zero, then (trial + 2*i) is composite.
364 ** Each prime used in the sieve costs one division of trial, and eliminates
365 ** one or more values from the search space. (3 eliminates 1/3 of the values
366 ** alone!) Each value left in the search space costs 1 or more modular
367 ** exponentations. So, these divisions are a bargain!
368 */
369 mp_err
mpp_sieve(mp_int * trial,const mp_digit * primes,mp_size nPrimes,unsigned char * sieve,mp_size nSieve)370 mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes,
371 unsigned char *sieve, mp_size nSieve)
372 {
373 mp_err res;
374 mp_digit rem;
375 mp_size ix;
376 unsigned long offset;
377
378 memset(sieve, 0, nSieve);
379
380 for (ix = 0; ix < nPrimes; ix++) {
381 mp_digit prime = primes[ix];
382 mp_size i;
383 if ((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY)
384 return res;
385
386 if (rem == 0) {
387 offset = 0;
388 } else {
389 offset = prime - rem;
390 }
391
392 for (i = offset; i < nSieve * 2; i += prime) {
393 if (i % 2 == 0) {
394 sieve[i / 2] = 1;
395 }
396 }
397 }
398
399 return MP_OKAY;
400 }
401
402 #define SIEVE_SIZE 32 * 1024
403
404 mp_err
mpp_make_prime(mp_int * start,mp_size nBits,mp_size strong,unsigned long * nTries)405 mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong,
406 unsigned long *nTries)
407 {
408 mp_digit np;
409 mp_err res;
410 unsigned int i = 0;
411 mp_int trial;
412 mp_int q;
413 mp_size num_tests;
414 unsigned char *sieve;
415
416 ARGCHK(start != 0, MP_BADARG);
417 ARGCHK(nBits > 16, MP_RANGE);
418
419 sieve = malloc(SIEVE_SIZE);
420 ARGCHK(sieve != NULL, MP_MEM);
421
422 MP_DIGITS(&trial) = 0;
423 MP_DIGITS(&q) = 0;
424 MP_CHECKOK(mp_init(&trial));
425 MP_CHECKOK(mp_init(&q));
426 /* values originally taken from table 4.4,
427 * HandBook of Applied Cryptography, augmented by FIPS-186
428 * requirements, Table C.2 and C.3 */
429 if (nBits >= 2000) {
430 num_tests = 3;
431 } else if (nBits >= 1536) {
432 num_tests = 4;
433 } else if (nBits >= 1024) {
434 num_tests = 5;
435 } else if (nBits >= 550) {
436 num_tests = 6;
437 } else if (nBits >= 450) {
438 num_tests = 7;
439 } else if (nBits >= 400) {
440 num_tests = 8;
441 } else if (nBits >= 350) {
442 num_tests = 9;
443 } else if (nBits >= 300) {
444 num_tests = 10;
445 } else if (nBits >= 250) {
446 num_tests = 20;
447 } else if (nBits >= 200) {
448 num_tests = 41;
449 } else if (nBits >= 100) {
450 num_tests = 38; /* funny anomaly in the FIPS tables, for aux primes, the
451 * required more iterations for larger aux primes */
452 } else
453 num_tests = 50;
454
455 if (strong)
456 --nBits;
457 MP_CHECKOK(mpl_set_bit(start, nBits - 1, 1));
458 MP_CHECKOK(mpl_set_bit(start, 0, 1));
459 for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) {
460 MP_CHECKOK(mpl_set_bit(start, i, 0));
461 }
462 /* start sieveing with prime value of 3. */
463 MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1,
464 sieve, SIEVE_SIZE));
465
466 #ifdef DEBUG_SIEVE
467 res = 0;
468 for (i = 0; i < SIEVE_SIZE; ++i) {
469 if (!sieve[i])
470 ++res;
471 }
472 fprintf(stderr, "sieve found %d potential primes.\n", res);
473 #define FPUTC(x, y) fputc(x, y)
474 #else
475 #define FPUTC(x, y)
476 #endif
477
478 res = MP_NO;
479 for (i = 0; i < SIEVE_SIZE; ++i) {
480 if (sieve[i]) /* this number is composite */
481 continue;
482 MP_CHECKOK(mp_add_d(start, 2 * i, &trial));
483 FPUTC('.', stderr);
484 /* run a Fermat test */
485 res = mpp_fermat(&trial, 2);
486 if (res != MP_OKAY) {
487 if (res == MP_NO)
488 continue; /* was composite */
489 goto CLEANUP;
490 }
491
492 FPUTC('+', stderr);
493 /* If that passed, run some Miller-Rabin tests */
494 res = mpp_pprime(&trial, num_tests);
495 if (res != MP_OKAY) {
496 if (res == MP_NO)
497 continue; /* was composite */
498 goto CLEANUP;
499 }
500 FPUTC('!', stderr);
501
502 if (!strong)
503 break; /* success !! */
504
505 /* At this point, we have strong evidence that our candidate
506 is itself prime. If we want a strong prime, we need now
507 to test q = 2p + 1 for primality...
508 */
509 MP_CHECKOK(mp_mul_2(&trial, &q));
510 MP_CHECKOK(mp_add_d(&q, 1, &q));
511
512 /* Test q for small prime divisors ... */
513 np = prime_tab_size;
514 res = mpp_divis_primes(&q, &np);
515 if (res == MP_YES) { /* is composite */
516 mp_clear(&q);
517 continue;
518 }
519 if (res != MP_NO)
520 goto CLEANUP;
521
522 /* And test with Fermat, as with its parent ... */
523 res = mpp_fermat(&q, 2);
524 if (res != MP_YES) {
525 mp_clear(&q);
526 if (res == MP_NO)
527 continue; /* was composite */
528 goto CLEANUP;
529 }
530
531 /* And test with Miller-Rabin, as with its parent ... */
532 res = mpp_pprime(&q, num_tests);
533 if (res != MP_YES) {
534 mp_clear(&q);
535 if (res == MP_NO)
536 continue; /* was composite */
537 goto CLEANUP;
538 }
539
540 /* If it passed, we've got a winner */
541 mp_exch(&q, &trial);
542 mp_clear(&q);
543 break;
544
545 } /* end of loop through sieved values */
546 if (res == MP_YES)
547 mp_exch(&trial, start);
548 CLEANUP:
549 mp_clear(&trial);
550 mp_clear(&q);
551 if (nTries)
552 *nTries += i;
553 if (sieve != NULL) {
554 memset(sieve, 0, SIEVE_SIZE);
555 free(sieve);
556 }
557 return res;
558 }
559
560 /*========================================================================*/
561 /*------------------------------------------------------------------------*/
562 /* Static functions visible only to the library internally */
563
564 /* {{{ s_mpp_divp(a, vec, size, which) */
565
566 /*
567 Test for divisibility by members of a vector of digits. Returns
568 MP_NO if a is not divisible by any of them; returns MP_YES and sets
569 'which' to the index of the offender, if it is. Will stop on the
570 first digit against which a is divisible.
571 */
572
573 mp_err
s_mpp_divp(mp_int * a,const mp_digit * vec,int size,int * which)574 s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which)
575 {
576 mp_err res;
577 mp_digit rem;
578
579 int ix;
580
581 for (ix = 0; ix < size; ix++) {
582 if ((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY)
583 return res;
584
585 if (rem == 0) {
586 if (which)
587 *which = ix;
588 return MP_YES;
589 }
590 }
591
592 return MP_NO;
593
594 } /* end s_mpp_divp() */
595
596 /* }}} */
597
598 /*------------------------------------------------------------------------*/
599 /* HERE THERE BE DRAGONS */
600