DXR is a code search and navigation tool aimed at making sense of large projects. It supports full-text and regex searches as well as structural queries.

Header

Untracked file

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622
/*
 *  mpprime.c
 *
 *  Utilities for finding and working with prime and pseudo-prime
 *  integers
 *
 * The contents of this file are subject to the Mozilla Public
 * License Version 1.1 (the "License"); you may not use this file
 * except in compliance with the License. You may obtain a copy of
 * the License at http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS
 * IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
 * implied. See the License for the specific language governing
 * rights and limitations under the License.
 *
 * The Original Code is the MPI Arbitrary Precision Integer Arithmetic
 * library.
 *
 * The Initial Developer of the Original Code is Michael J. Fromberger.
 * Portions created by Michael J. Fromberger are 
 * Copyright (C) 1997, 1998, 1999, 2000 Michael J. Fromberger. 
 * All Rights Reserved.
 *
 * Contributor(s):
 *	Netscape Communications Corporation
 *
 * Alternatively, the contents of this file may be used under the
 * terms of the GNU General Public License Version 2 or later (the
 * "GPL"), in which case the provisions of the GPL are applicable
 * instead of those above.  If you wish to allow use of your
 * version of this file only under the terms of the GPL and not to
 * allow others to use your version of this file under the MPL,
 * indicate your decision by deleting the provisions above and
 * replace them with the notice and other provisions required by
 * the GPL.  If you do not delete the provisions above, a recipient
 * may use your version of this file under either the MPL or the GPL.
 */

#include "mpi-priv.h"
#include "mpprime.h"
#include "mplogic.h"
#include <stdlib.h>
#include <string.h>

#define SMALL_TABLE 0 /* determines size of hard-wired prime table */

#define RANDOM() rand()

#include "primes.c"  /* pull in the prime digit table */

/* 
   Test if any of a given vector of digits divides a.  If not, MP_NO
   is returned; otherwise, MP_YES is returned and 'which' is set to
   the index of the integer in the vector which divided a.
 */
mp_err    s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which);

/* {{{ mpp_divis(a, b) */

/*
  mpp_divis(a, b)

  Returns MP_YES if a is divisible by b, or MP_NO if it is not.
 */

mp_err  mpp_divis(mp_int *a, mp_int *b)
{
  mp_err  res;
  mp_int  rem;

  if((res = mp_init(&rem)) != MP_OKAY)
    return res;

  if((res = mp_mod(a, b, &rem)) != MP_OKAY)
    goto CLEANUP;

  if(mp_cmp_z(&rem) == 0)
    res = MP_YES;
  else
    res = MP_NO;

CLEANUP:
  mp_clear(&rem);
  return res;

} /* end mpp_divis() */

/* }}} */

/* {{{ mpp_divis_d(a, d) */

/*
  mpp_divis_d(a, d)

  Return MP_YES if a is divisible by d, or MP_NO if it is not.
 */

mp_err  mpp_divis_d(mp_int *a, mp_digit d)
{
  mp_err     res;
  mp_digit   rem;

  ARGCHK(a != NULL, MP_BADARG);

  if(d == 0)
    return MP_NO;

  if((res = mp_mod_d(a, d, &rem)) != MP_OKAY)
    return res;

  if(rem == 0)
    return MP_YES;
  else
    return MP_NO;

} /* end mpp_divis_d() */

/* }}} */

/* {{{ mpp_random(a) */

/*
  mpp_random(a)

  Assigns a random value to a.  This value is generated using the
  standard C library's rand() function, so it should not be used for
  cryptographic purposes, but it should be fine for primality testing,
  since all we really care about there is good statistical properties.

  As many digits as a currently has are filled with random digits.
 */

mp_err  mpp_random(mp_int *a)

{
  mp_digit  next = 0;
  unsigned int       ix, jx;

  ARGCHK(a != NULL, MP_BADARG);

  for(ix = 0; ix < USED(a); ix++) {
    for(jx = 0; jx < sizeof(mp_digit); jx++) {
      next = (next << CHAR_BIT) | (RANDOM() & UCHAR_MAX);
    }
    DIGIT(a, ix) = next;
  }

  return MP_OKAY;

} /* end mpp_random() */

/* }}} */

/* {{{ mpp_random_size(a, prec) */

mp_err  mpp_random_size(mp_int *a, mp_size prec)
{
  mp_err   res;

  ARGCHK(a != NULL && prec > 0, MP_BADARG);
  
  if((res = s_mp_pad(a, prec)) != MP_OKAY)
    return res;

  return mpp_random(a);

} /* end mpp_random_size() */

/* }}} */

/* {{{ mpp_divis_vector(a, vec, size, which) */

/*
  mpp_divis_vector(a, vec, size, which)

  Determines if a is divisible by any of the 'size' digits in vec.
  Returns MP_YES and sets 'which' to the index of the offending digit,
  if it is; returns MP_NO if it is not.
 */

mp_err  mpp_divis_vector(mp_int *a, const mp_digit *vec, int size, int *which)
{
  ARGCHK(a != NULL && vec != NULL && size > 0, MP_BADARG);
  
  return s_mpp_divp(a, vec, size, which);

} /* end mpp_divis_vector() */

/* }}} */

/* {{{ mpp_divis_primes(a, np) */

/*
  mpp_divis_primes(a, np)

  Test whether a is divisible by any of the first 'np' primes.  If it
  is, returns MP_YES and sets *np to the value of the digit that did
  it.  If not, returns MP_NO.
 */
mp_err  mpp_divis_primes(mp_int *a, mp_digit *np)
{
  int     size, which;
  mp_err  res;

  ARGCHK(a != NULL && np != NULL, MP_BADARG);

  size = (int)*np;
  if(size > prime_tab_size)
    size = prime_tab_size;

  res = mpp_divis_vector(a, prime_tab, size, &which);
  if(res == MP_YES) 
    *np = prime_tab[which];

  return res;

} /* end mpp_divis_primes() */

/* }}} */

/* {{{ mpp_fermat(a, w) */

/*
  Using w as a witness, try pseudo-primality testing based on Fermat's
  little theorem.  If a is prime, and (w, a) = 1, then w^a == w (mod
  a).  So, we compute z = w^a (mod a) and compare z to w; if they are
  equal, the test passes and we return MP_YES.  Otherwise, we return
  MP_NO.
 */
mp_err  mpp_fermat(mp_int *a, mp_digit w)
{
  mp_int  base, test;
  mp_err  res;
  
  if((res = mp_init(&base)) != MP_OKAY)
    return res;

  mp_set(&base, w);

  if((res = mp_init(&test)) != MP_OKAY)
    goto TEST;

  /* Compute test = base^a (mod a) */
  if((res = mp_exptmod(&base, a, a, &test)) != MP_OKAY)
    goto CLEANUP;

  
  if(mp_cmp(&base, &test) == 0)
    res = MP_YES;
  else
    res = MP_NO;

 CLEANUP:
  mp_clear(&test);
 TEST:
  mp_clear(&base);

  return res;

} /* end mpp_fermat() */

/* }}} */

/*
  Perform the fermat test on each of the primes in a list until
  a) one of them shows a is not prime, or 
  b) the list is exhausted.
  Returns:  MP_YES if it passes tests.
	    MP_NO  if fermat test reveals it is composite
	    Some MP error code if some other error occurs.
 */
mp_err mpp_fermat_list(mp_int *a, const mp_digit *primes, mp_size nPrimes)
{
  mp_err rv = MP_YES;

  while (nPrimes-- > 0 && rv == MP_YES) {
    rv = mpp_fermat(a, *primes++);
  }
  return rv;
}

/* {{{ mpp_pprime(a, nt) */

/*
  mpp_pprime(a, nt)

  Performs nt iteration of the Miller-Rabin probabilistic primality
  test on a.  Returns MP_YES if the tests pass, MP_NO if one fails.
  If MP_NO is returned, the number is definitely composite.  If MP_YES
  is returned, it is probably prime (but that is not guaranteed).
 */

mp_err  mpp_pprime(mp_int *a, int nt)
{
  mp_err   res;
  mp_int   x, amo, m, z;	/* "amo" = "a minus one" */
  int      iter;
  unsigned int jx;
  mp_size  b;

  ARGCHK(a != NULL, MP_BADARG);

  MP_DIGITS(&x) = 0;
  MP_DIGITS(&amo) = 0;
  MP_DIGITS(&m) = 0;
  MP_DIGITS(&z) = 0;

  /* Initialize temporaries... */
  MP_CHECKOK( mp_init(&amo));
  /* Compute amo = a - 1 for what follows...    */
  MP_CHECKOK( mp_sub_d(a, 1, &amo) );

  b = mp_trailing_zeros(&amo);
  if (!b) { /* a was even ? */
    res = MP_NO;
    goto CLEANUP;
  }

  MP_CHECKOK( mp_init_size(&x, MP_USED(a)) );
  MP_CHECKOK( mp_init(&z) );
  MP_CHECKOK( mp_init(&m) );
  MP_CHECKOK( mp_div_2d(&amo, b, &m, 0) );

  /* Do the test nt times... */
  for(iter = 0; iter < nt; iter++) {

    /* Choose a random value for x < a          */
    s_mp_pad(&x, USED(a));
    mpp_random(&x);
    MP_CHECKOK( mp_mod(&x, a, &x) );

    /* Compute z = (x ** m) mod a               */
    MP_CHECKOK( mp_exptmod(&x, &m, a, &z) );
    
    if(mp_cmp_d(&z, 1) == 0 || mp_cmp(&z, &amo) == 0) {
      res = MP_YES;
      continue;
    }
    
    res = MP_NO;  /* just in case the following for loop never executes. */
    for (jx = 1; jx < b; jx++) {
      /* z = z^2 (mod a) */
      MP_CHECKOK( mp_sqrmod(&z, a, &z) );
      res = MP_NO;	/* previous line set res to MP_YES */

      if(mp_cmp_d(&z, 1) == 0) {
	break;
      }
      if(mp_cmp(&z, &amo) == 0) {
	res = MP_YES;
	break;
      } 
    } /* end testing loop */

    /* If the test passes, we will continue iterating, but a failed
       test means the candidate is definitely NOT prime, so we will
       immediately break out of this loop
     */
    if(res == MP_NO)
      break;

  } /* end iterations loop */
  
CLEANUP:
  mp_clear(&m);
  mp_clear(&z);
  mp_clear(&x);
  mp_clear(&amo);
  return res;

} /* end mpp_pprime() */

/* }}} */

/* Produce table of composites from list of primes and trial value.  
** trial must be odd. List of primes must not include 2.
** sieve should have dimension >= MAXPRIME/2, where MAXPRIME is largest 
** prime in list of primes.  After this function is finished,
** if sieve[i] is non-zero, then (trial + 2*i) is composite.
** Each prime used in the sieve costs one division of trial, and eliminates
** one or more values from the search space. (3 eliminates 1/3 of the values
** alone!)  Each value left in the search space costs 1 or more modular 
** exponentations.  So, these divisions are a bargain!
*/
mp_err mpp_sieve(mp_int *trial, const mp_digit *primes, mp_size nPrimes, 
		 unsigned char *sieve, mp_size nSieve)
{
  mp_err       res;
  mp_digit     rem;
  mp_size      ix;
  unsigned long offset;

  memset(sieve, 0, nSieve);

  for(ix = 0; ix < nPrimes; ix++) {
    mp_digit prime = primes[ix];
    mp_size  i;
    if((res = mp_mod_d(trial, prime, &rem)) != MP_OKAY) 
      return res;

    if (rem == 0) {
      offset = 0;
    } else {
      offset = prime - (rem / 2);
    }
    for (i = offset; i < nSieve ; i += prime) {
      sieve[i] = 1;
    }
  }

  return MP_OKAY;
}

#define SIEVE_SIZE 32*1024

mp_err mpp_make_prime(mp_int *start, mp_size nBits, mp_size strong,
		      unsigned long * nTries)
{
  mp_digit      np;
  mp_err        res;
  int           i	= 0;
  mp_int        trial;
  mp_int        q;
  mp_size       num_tests;
  /*
   * Always make sieve the last variabale allocated so that 
   * Mac builds don't break by adding an extra variable
   * on the stack. -javi
   */
#if defined(macintosh) || defined (XP_OS2) \
    || (defined(HPUX) && defined(__ia64))
  unsigned char *sieve;
  
  sieve = malloc(SIEVE_SIZE);
  ARGCHK(sieve != NULL, MP_MEM);
#else
  unsigned char sieve[SIEVE_SIZE];
#endif  

  ARGCHK(start != 0, MP_BADARG);
  ARGCHK(nBits > 16, MP_RANGE);

  MP_DIGITS(&trial) = 0;
  MP_DIGITS(&q) = 0;
  MP_CHECKOK( mp_init(&trial) );
  MP_CHECKOK( mp_init(&q)     );
  /* values taken from table 4.4, HandBook of Applied Cryptography */
  if (nBits >= 1300) {
    num_tests = 2;
  } else if (nBits >= 850) {
    num_tests = 3;
  } else if (nBits >= 650) {
    num_tests = 4;
  } else if (nBits >= 550) {
    num_tests = 5;
  } else if (nBits >= 450) {
    num_tests = 6;
  } else if (nBits >= 400) {
    num_tests = 7;
  } else if (nBits >= 350) {
    num_tests = 8;
  } else if (nBits >= 300) {
    num_tests = 9;
  } else if (nBits >= 250) {
    num_tests = 12;
  } else if (nBits >= 200) {
    num_tests = 15;
  } else if (nBits >= 150) {
    num_tests = 18;
  } else if (nBits >= 100) {
    num_tests = 27;
  } else
    num_tests = 50;

  if (strong) 
    --nBits;
  MP_CHECKOK( mpl_set_bit(start, nBits - 1, 1) );
  MP_CHECKOK( mpl_set_bit(start,         0, 1) );
  for (i = mpl_significant_bits(start) - 1; i >= nBits; --i) {
    MP_CHECKOK( mpl_set_bit(start, i, 0) );
  }
  /* start sieveing with prime value of 3. */
  MP_CHECKOK(mpp_sieve(start, prime_tab + 1, prime_tab_size - 1, 
		       sieve, SIEVE_SIZE) );

#ifdef DEBUG_SIEVE
  res = 0;
  for (i = 0; i < SIEVE_SIZE; ++i) {
    if (!sieve[i])
      ++res;
  }
  fprintf(stderr,"sieve found %d potential primes.\n", res);
#define FPUTC(x,y) fputc(x,y)
#else
#define FPUTC(x,y) 
#endif

  res = MP_NO;
  for(i = 0; i < SIEVE_SIZE; ++i) {
    if (sieve[i])	/* this number is composite */
      continue;
    MP_CHECKOK( mp_add_d(start, 2 * i, &trial) );
    FPUTC('.', stderr);
    /* run a Fermat test */
    res = mpp_fermat(&trial, 2);
    if (res != MP_OKAY) {
      if (res == MP_NO)
	continue;	/* was composite */
      goto CLEANUP;
    }
      
    FPUTC('+', stderr);
    /* If that passed, run some Miller-Rabin tests	*/
    res = mpp_pprime(&trial, num_tests);
    if (res != MP_OKAY) {
      if (res == MP_NO)
	continue;	/* was composite */
      goto CLEANUP;
    }
    FPUTC('!', stderr);

    if (!strong) 
      break;	/* success !! */

    /* At this point, we have strong evidence that our candidate
       is itself prime.  If we want a strong prime, we need now
       to test q = 2p + 1 for primality...
     */
    MP_CHECKOK( mp_mul_2(&trial, &q) );
    MP_CHECKOK( mp_add_d(&q, 1, &q)  );

    /* Test q for small prime divisors ... */
    np = prime_tab_size;
    res = mpp_divis_primes(&q, &np);
    if (res == MP_YES) { /* is composite */
      mp_clear(&q);
      continue;
    }
    if (res != MP_NO) 
      goto CLEANUP;

    /* And test with Fermat, as with its parent ... */
    res = mpp_fermat(&q, 2);
    if (res != MP_YES) {
      mp_clear(&q);
      if (res == MP_NO)
	continue;	/* was composite */
      goto CLEANUP;
    }

    /* And test with Miller-Rabin, as with its parent ... */
    res = mpp_pprime(&q, num_tests);
    if (res != MP_YES) {
      mp_clear(&q);
      if (res == MP_NO)
	continue;	/* was composite */
      goto CLEANUP;
    }

    /* If it passed, we've got a winner */
    mp_exch(&q, &trial);
    mp_clear(&q);
    break;

  } /* end of loop through sieved values */
  if (res == MP_YES) 
    mp_exch(&trial, start);
CLEANUP:
  mp_clear(&trial);
  mp_clear(&q);
  if (nTries)
    *nTries += i;
#if defined(macintosh) || defined(XP_OS2) \
    || (defined(HPUX) && defined(__ia64))
  if (sieve != NULL) {
  	memset(sieve, 0, SIEVE_SIZE);
  	free (sieve);
  }
#endif    
  return res;
}

/*========================================================================*/
/*------------------------------------------------------------------------*/
/* Static functions visible only to the library internally                */

/* {{{ s_mpp_divp(a, vec, size, which) */

/* 
   Test for divisibility by members of a vector of digits.  Returns
   MP_NO if a is not divisible by any of them; returns MP_YES and sets
   'which' to the index of the offender, if it is.  Will stop on the
   first digit against which a is divisible.
 */

mp_err    s_mpp_divp(mp_int *a, const mp_digit *vec, int size, int *which)
{
  mp_err    res;
  mp_digit  rem;

  int     ix;

  for(ix = 0; ix < size; ix++) {
    if((res = mp_mod_d(a, vec[ix], &rem)) != MP_OKAY) 
      return res;

    if(rem == 0) {
      if(which)
	*which = ix;
      return MP_YES;
    }
  }

  return MP_NO;

} /* end s_mpp_divp() */

/* }}} */

/*------------------------------------------------------------------------*/
/* HERE THERE BE DRAGONS                                                  */