The Tor Browser / file comparison

comparison: security/nss/lib/freebl/mpi/mpprime.c

security/nss/lib/freebl/mpi/mpprime.c

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

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial