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