version 1.1.1.1, 2000/01/10 15:35:27 |
version 1.1.1.4, 2003/08/25 16:06:33 |
|
|
An implementation of the probabilistic primality test found in Knuth's |
An implementation of the probabilistic primality test found in Knuth's |
Seminumerical Algorithms book. If the function mpz_probab_prime_p() |
Seminumerical Algorithms book. If the function mpz_probab_prime_p() |
returns 0 then n is not prime. If it returns 1, then n is 'probably' |
returns 0 then n is not prime. If it returns 1, then n is 'probably' |
prime. The probability of a false positive is (1/4)**reps, where |
prime. If it returns 2, n is surely prime. The probability of a false |
reps is the number of internal passes of the probabilistic algorithm. |
positive is (1/4)**reps, where reps is the number of internal passes of the |
Knuth indicates that 25 passes are reasonable. |
probabilistic algorithm. Knuth indicates that 25 passes are reasonable. |
|
|
Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc. |
Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software |
Contributed by John Amanatides. |
Foundation, Inc. Miller-Rabin code contributed by John Amanatides. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
The GNU MP Library is free software; you can redistribute it and/or modify |
The GNU MP Library is free software; you can redistribute it and/or modify |
it under the terms of the GNU Library General Public License as published by |
it under the terms of the GNU Lesser General Public License as published by |
the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
option) any later version. |
option) any later version. |
|
|
The GNU MP Library is distributed in the hope that it will be useful, but |
The GNU MP Library is distributed in the hope that it will be useful, but |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
License for more details. |
License for more details. |
|
|
You should have received a copy of the GNU Library General Public License |
You should have received a copy of the GNU Lesser General Public License |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
|
|
#include "gmp.h" |
#include "gmp.h" |
|
#include "gmp-impl.h" |
|
#include "longlong.h" |
|
|
static int |
static int isprime _PROTO ((unsigned long int t)); |
possibly_prime (n, n_minus_1, x, y, q, k) |
|
mpz_srcptr n; |
int |
mpz_srcptr n_minus_1; |
mpz_probab_prime_p (mpz_srcptr n, int reps) |
mpz_ptr x; |
|
mpz_ptr y; |
|
mpz_srcptr q; |
|
unsigned long int k; |
|
{ |
{ |
unsigned long int i; |
mp_limb_t r; |
|
|
/* find random x s.t. 1 < x < n */ |
/* Handle small and negative n. */ |
do |
if (mpz_cmp_ui (n, 1000000L) <= 0) |
{ |
{ |
mpz_random (x, mpz_size (n)); |
int is_prime; |
mpz_mmod (x, x, n); |
if (mpz_sgn (n) < 0) |
|
{ |
|
/* Negative number. Negate and call ourselves. */ |
|
mpz_t n2; |
|
mpz_init (n2); |
|
mpz_neg (n2, n); |
|
is_prime = mpz_probab_prime_p (n2, reps); |
|
mpz_clear (n2); |
|
return is_prime; |
|
} |
|
is_prime = isprime (mpz_get_ui (n)); |
|
return is_prime ? 2 : 0; |
} |
} |
while (mpz_cmp_ui (x, 1L) <= 0); |
|
|
|
mpz_powm (y, x, q, n); |
/* If n is now even, it is not a prime. */ |
|
if ((mpz_get_ui (n) & 1) == 0) |
|
return 0; |
|
|
if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0) |
#if defined (PP) |
return 1; |
/* Check if n has small factors. */ |
|
#if defined (PP_INVERTED) |
for (i = 1; i < k; i++) |
r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP, |
{ |
(mp_limb_t) PP_INVERTED); |
mpz_powm_ui (y, y, 2L, n); |
|
if (mpz_cmp (y, n_minus_1) == 0) |
|
return 1; |
|
if (mpz_cmp_ui (y, 1L) == 0) |
|
return 0; |
|
} |
|
return 0; |
|
} |
|
|
|
int |
|
#if __STDC__ |
|
mpz_probab_prime_p (mpz_srcptr m, int reps) |
|
#else |
#else |
mpz_probab_prime_p (m, reps) |
r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); |
mpz_srcptr m; |
|
int reps; |
|
#endif |
#endif |
{ |
if (r % 3 == 0 |
mpz_t n, n_minus_1, x, y, q; |
#if BITS_PER_MP_LIMB >= 4 |
int i, is_prime; |
|| r % 5 == 0 |
unsigned long int k; |
#endif |
|
#if BITS_PER_MP_LIMB >= 8 |
mpz_init (n); |
|| r % 7 == 0 |
/* Take the absolute value of M, to handle positive and negative primes. */ |
#endif |
mpz_abs (n, m); |
#if BITS_PER_MP_LIMB >= 16 |
|
|| r % 11 == 0 || r % 13 == 0 |
if (mpz_cmp_ui (n, 3L) <= 0) |
#endif |
|
#if BITS_PER_MP_LIMB >= 32 |
|
|| r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 |
|
#endif |
|
#if BITS_PER_MP_LIMB >= 64 |
|
|| r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0 |
|
|| r % 47 == 0 || r % 53 == 0 |
|
#endif |
|
) |
{ |
{ |
mpz_clear (n); |
return 0; |
return mpz_cmp_ui (n, 1L) > 0; |
|
} |
} |
|
#endif /* PP */ |
|
|
if ((mpz_get_ui (n) & 1) == 0) |
/* Do more dividing. We collect small primes, using umul_ppmm, until we |
{ |
overflow a single limb. We divide our number by the small primes product, |
mpz_clear (n); |
and look for factors in the remainder. */ |
return 0; /* even */ |
{ |
} |
unsigned long int ln2; |
|
unsigned long int q; |
|
mp_limb_t p1, p0, p; |
|
unsigned int primes[15]; |
|
int nprimes; |
|
|
mpz_init (n_minus_1); |
nprimes = 0; |
mpz_sub_ui (n_minus_1, n, 1L); |
p = 1; |
mpz_init (x); |
ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; |
mpz_init (y); |
for (q = PP_FIRST_OMITTED; q < ln2; q += 2) |
|
{ |
|
if (isprime (q)) |
|
{ |
|
umul_ppmm (p1, p0, p, q); |
|
if (p1 != 0) |
|
{ |
|
r = mpn_mod_1 (PTR(n), SIZ(n), p); |
|
while (--nprimes >= 0) |
|
if (r % primes[nprimes] == 0) |
|
{ |
|
ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0); |
|
return 0; |
|
} |
|
p = q; |
|
nprimes = 0; |
|
} |
|
else |
|
{ |
|
p = p0; |
|
} |
|
primes[nprimes++] = q; |
|
} |
|
} |
|
} |
|
|
/* find q and k, s.t. n = 1 + 2**k * q */ |
/* Perform a number of Miller-Rabin tests. */ |
mpz_init_set (q, n_minus_1); |
return mpz_millerrabin (n, reps); |
k = mpz_scan1 (q, 0); |
} |
mpz_tdiv_q_2exp (q, q, k); |
|
|
|
is_prime = 1; |
static int |
for (i = 0; i < reps && is_prime; i++) |
isprime (unsigned long int t) |
is_prime &= possibly_prime (n, n_minus_1, x, y, q, k); |
{ |
|
unsigned long int q, r, d; |
|
|
mpz_clear (n_minus_1); |
if (t < 3 || (t & 1) == 0) |
mpz_clear (n); |
return t == 2; |
mpz_clear (x); |
|
mpz_clear (y); |
for (d = 3, r = 1; r != 0; d += 2) |
mpz_clear (q); |
{ |
return is_prime; |
q = t / d; |
|
r = t - q * d; |
|
if (q < d) |
|
return 1; |
|
} |
|
return 0; |
} |
} |