Annotation of OpenXM/src/kan96xx/gmp-2.0.2/mpz/pprime_p.c, Revision 1.1
1.1 ! maekawa 1: /* mpz_probab_prime_p --
! 2: An implementation of the probabilistic primality test found in Knuth's
! 3: Seminumerical Algorithms book. If the function mpz_probab_prime_p()
! 4: returns 0 then n is not prime. If it returns 1, then n is 'probably'
! 5: prime. The probability of a false positive is (1/4)**reps, where
! 6: reps is the number of internal passes of the probabilistic algorithm.
! 7: Knuth indicates that 25 passes are reasonable.
! 8:
! 9: Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc.
! 10: Contributed by John Amanatides.
! 11:
! 12: This file is part of the GNU MP Library.
! 13:
! 14: The GNU MP Library is free software; you can redistribute it and/or modify
! 15: it under the terms of the GNU Library General Public License as published by
! 16: the Free Software Foundation; either version 2 of the License, or (at your
! 17: option) any later version.
! 18:
! 19: The GNU MP Library is distributed in the hope that it will be useful, but
! 20: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 21: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 22: License for more details.
! 23:
! 24: You should have received a copy of the GNU Library General Public License
! 25: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 26: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 27: MA 02111-1307, USA. */
! 28:
! 29: #include "gmp.h"
! 30:
! 31: static int
! 32: possibly_prime (n, n_minus_1, x, y, q, k)
! 33: mpz_srcptr n;
! 34: mpz_srcptr n_minus_1;
! 35: mpz_ptr x;
! 36: mpz_ptr y;
! 37: mpz_srcptr q;
! 38: unsigned long int k;
! 39: {
! 40: unsigned long int i;
! 41:
! 42: /* find random x s.t. 1 < x < n */
! 43: do
! 44: {
! 45: mpz_random (x, mpz_size (n));
! 46: mpz_mmod (x, x, n);
! 47: }
! 48: while (mpz_cmp_ui (x, 1L) <= 0);
! 49:
! 50: mpz_powm (y, x, q, n);
! 51:
! 52: if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0)
! 53: return 1;
! 54:
! 55: for (i = 1; i < k; i++)
! 56: {
! 57: mpz_powm_ui (y, y, 2L, n);
! 58: if (mpz_cmp (y, n_minus_1) == 0)
! 59: return 1;
! 60: if (mpz_cmp_ui (y, 1L) == 0)
! 61: return 0;
! 62: }
! 63: return 0;
! 64: }
! 65:
! 66: int
! 67: #if __STDC__
! 68: mpz_probab_prime_p (mpz_srcptr m, int reps)
! 69: #else
! 70: mpz_probab_prime_p (m, reps)
! 71: mpz_srcptr m;
! 72: int reps;
! 73: #endif
! 74: {
! 75: mpz_t n, n_minus_1, x, y, q;
! 76: int i, is_prime;
! 77: unsigned long int k;
! 78:
! 79: mpz_init (n);
! 80: /* Take the absolute value of M, to handle positive and negative primes. */
! 81: mpz_abs (n, m);
! 82:
! 83: if (mpz_cmp_ui (n, 3L) <= 0)
! 84: {
! 85: mpz_clear (n);
! 86: return mpz_cmp_ui (n, 1L) > 0;
! 87: }
! 88:
! 89: if ((mpz_get_ui (n) & 1) == 0)
! 90: {
! 91: mpz_clear (n);
! 92: return 0; /* even */
! 93: }
! 94:
! 95: mpz_init (n_minus_1);
! 96: mpz_sub_ui (n_minus_1, n, 1L);
! 97: mpz_init (x);
! 98: mpz_init (y);
! 99:
! 100: /* find q and k, s.t. n = 1 + 2**k * q */
! 101: mpz_init_set (q, n_minus_1);
! 102: k = mpz_scan1 (q, 0);
! 103: mpz_tdiv_q_2exp (q, q, k);
! 104:
! 105: is_prime = 1;
! 106: for (i = 0; i < reps && is_prime; i++)
! 107: is_prime &= possibly_prime (n, n_minus_1, x, y, q, k);
! 108:
! 109: mpz_clear (n_minus_1);
! 110: mpz_clear (n);
! 111: mpz_clear (x);
! 112: mpz_clear (y);
! 113: mpz_clear (q);
! 114: return is_prime;
! 115: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>