version 1.1, 2000/01/10 15:35:27 |
version 1.1.1.3, 2000/12/01 05:45:11 |
|
|
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 (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 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) |
static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps)); |
|
|
|
int |
|
#if __STDC__ |
|
mpz_probab_prime_p (mpz_srcptr n, int reps) |
|
#else |
|
mpz_probab_prime_p (n, reps) |
mpz_srcptr n; |
mpz_srcptr n; |
mpz_srcptr n_minus_1; |
int reps; |
mpz_ptr x; |
#endif |
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) |
/* Check if n has small factors. */ |
return 1; |
if (UDIV_TIME > (2 * UMUL_TIME + 6)) |
|
r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED); |
|
else |
|
r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); |
|
if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0 |
|
|| r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 |
|
#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 |
|
) |
|
{ |
|
return 0; |
|
} |
|
|
for (i = 1; i < k; i++) |
/* 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, |
|
and look for factors in the remainder. */ |
|
{ |
|
unsigned long int ln2; |
|
unsigned long int q; |
|
mp_limb_t p1, p0, p; |
|
unsigned int primes[15]; |
|
int nprimes; |
|
|
|
nprimes = 0; |
|
p = 1; |
|
ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; |
|
for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; 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) |
|
{ |
|
if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0) |
|
abort (); |
|
return 0; |
|
} |
|
p = q; |
|
nprimes = 0; |
|
} |
|
else |
|
{ |
|
p = p0; |
|
} |
|
primes[nprimes++] = q; |
|
} |
|
} |
|
} |
|
|
|
/* Perform a number of Miller-Rabin tests. */ |
|
return mpz_millerrabin (n, reps); |
|
} |
|
|
|
static int |
|
#if __STDC__ |
|
isprime (unsigned long int t) |
|
#else |
|
isprime (t) |
|
unsigned long int t; |
|
#endif |
|
{ |
|
unsigned long int q, r, d; |
|
|
|
if (t < 3 || (t & 1) == 0) |
|
return t == 2; |
|
|
|
for (d = 3, r = 1; r != 0; d += 2) |
{ |
{ |
mpz_powm_ui (y, y, 2L, n); |
q = t / d; |
if (mpz_cmp (y, n_minus_1) == 0) |
r = t - q * d; |
|
if (q < d) |
return 1; |
return 1; |
if (mpz_cmp_ui (y, 1L) == 0) |
|
return 0; |
|
} |
} |
return 0; |
return 0; |
} |
} |
|
|
int |
static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1, |
|
mpz_ptr x, mpz_ptr y, |
|
mpz_srcptr q, unsigned long int k)); |
|
|
|
static int |
#if __STDC__ |
#if __STDC__ |
mpz_probab_prime_p (mpz_srcptr m, int reps) |
mpz_millerrabin (mpz_srcptr n, int reps) |
#else |
#else |
mpz_probab_prime_p (m, reps) |
mpz_millerrabin (n, reps) |
mpz_srcptr m; |
mpz_srcptr n; |
int reps; |
int reps; |
#endif |
#endif |
{ |
{ |
mpz_t n, n_minus_1, x, y, q; |
int r; |
int i, is_prime; |
mpz_t nm1, x, y, q; |
unsigned long int k; |
unsigned long int k; |
|
gmp_randstate_t rstate; |
|
int is_prime; |
|
TMP_DECL (marker); |
|
TMP_MARK (marker); |
|
|
mpz_init (n); |
MPZ_TMP_INIT (nm1, SIZ (n) + 1); |
/* Take the absolute value of M, to handle positive and negative primes. */ |
mpz_sub_ui (nm1, n, 1L); |
mpz_abs (n, m); |
|
|
|
if (mpz_cmp_ui (n, 3L) <= 0) |
MPZ_TMP_INIT (x, SIZ (n)); |
{ |
MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ |
mpz_clear (n); |
|
return mpz_cmp_ui (n, 1L) > 0; |
|
} |
|
|
|
if ((mpz_get_ui (n) & 1) == 0) |
/* Perform a Fermat test. */ |
|
mpz_set_ui (x, 210L); |
|
mpz_powm (y, x, nm1, n); |
|
if (mpz_cmp_ui (y, 1L) != 0) |
{ |
{ |
mpz_clear (n); |
TMP_FREE (marker); |
return 0; /* even */ |
return 0; |
} |
} |
|
|
mpz_init (n_minus_1); |
MPZ_TMP_INIT (q, SIZ (n)); |
mpz_sub_ui (n_minus_1, n, 1L); |
|
mpz_init (x); |
|
mpz_init (y); |
|
|
|
/* find q and k, s.t. n = 1 + 2**k * q */ |
/* Find q and k, where q is odd and n = 1 + 2**k * q. */ |
mpz_init_set (q, n_minus_1); |
k = mpz_scan1 (nm1, 0L); |
k = mpz_scan1 (q, 0); |
mpz_tdiv_q_2exp (q, nm1, k); |
mpz_tdiv_q_2exp (q, q, k); |
|
|
|
|
gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L); |
|
|
is_prime = 1; |
is_prime = 1; |
for (i = 0; i < reps && is_prime; i++) |
for (r = 0; r < reps && is_prime; r++) |
is_prime &= possibly_prime (n, n_minus_1, x, y, q, k); |
{ |
|
do |
|
mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1); |
|
while (mpz_cmp_ui (x, 1L) <= 0); |
|
|
mpz_clear (n_minus_1); |
is_prime = millerrabin (n, nm1, x, y, q, k); |
mpz_clear (n); |
} |
mpz_clear (x); |
|
mpz_clear (y); |
gmp_randclear (rstate); |
mpz_clear (q); |
|
|
TMP_FREE (marker); |
return is_prime; |
return is_prime; |
|
} |
|
|
|
static int |
|
#if __STDC__ |
|
millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y, |
|
mpz_srcptr q, unsigned long int k) |
|
#else |
|
millerrabin (n, nm1, x, y, q, k) |
|
mpz_srcptr n; |
|
mpz_srcptr nm1; |
|
mpz_ptr x; |
|
mpz_ptr y; |
|
mpz_srcptr q; |
|
unsigned long int k; |
|
#endif |
|
{ |
|
unsigned long int i; |
|
|
|
mpz_powm (y, x, q, n); |
|
|
|
if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0) |
|
return 1; |
|
|
|
for (i = 1; i < k; i++) |
|
{ |
|
mpz_powm_ui (y, y, 2L, n); |
|
if (mpz_cmp (y, nm1) == 0) |
|
return 1; |
|
if (mpz_cmp_ui (y, 1L) == 0) |
|
return 0; |
|
} |
|
return 0; |
} |
} |