version 1.1.1.3, 2000/12/01 05:45:11 |
version 1.1.1.4, 2003/08/25 16:06:33 |
|
|
positive is (1/4)**reps, where reps is the number of internal passes of the |
positive is (1/4)**reps, where reps is the number of internal passes of the |
probabilistic algorithm. Knuth indicates that 25 passes are reasonable. |
probabilistic algorithm. Knuth indicates that 25 passes are reasonable. |
|
|
Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software |
Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software |
Foundation, Inc. Miller-Rabin code 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. |
Line 31 MA 02111-1307, USA. */ |
|
Line 31 MA 02111-1307, USA. */ |
|
#include "longlong.h" |
#include "longlong.h" |
|
|
static int isprime _PROTO ((unsigned long int t)); |
static int isprime _PROTO ((unsigned long int t)); |
static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps)); |
|
|
|
int |
int |
#if __STDC__ |
|
mpz_probab_prime_p (mpz_srcptr n, int reps) |
mpz_probab_prime_p (mpz_srcptr n, int reps) |
#else |
|
mpz_probab_prime_p (n, reps) |
|
mpz_srcptr n; |
|
int reps; |
|
#endif |
|
{ |
{ |
mp_limb_t r; |
mp_limb_t r; |
|
|
Line 66 mpz_probab_prime_p (n, reps) |
|
Line 59 mpz_probab_prime_p (n, reps) |
|
if ((mpz_get_ui (n) & 1) == 0) |
if ((mpz_get_ui (n) & 1) == 0) |
return 0; |
return 0; |
|
|
|
#if defined (PP) |
/* Check if n has small factors. */ |
/* Check if n has small factors. */ |
if (UDIV_TIME > (2 * UMUL_TIME + 6)) |
#if defined (PP_INVERTED) |
r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED); |
r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP, |
else |
(mp_limb_t) PP_INVERTED); |
r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); |
#else |
if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0 |
r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP); |
|
#endif |
|
if (r % 3 == 0 |
|
#if BITS_PER_MP_LIMB >= 4 |
|
|| r % 5 == 0 |
|
#endif |
|
#if BITS_PER_MP_LIMB >= 8 |
|
|| r % 7 == 0 |
|
#endif |
|
#if BITS_PER_MP_LIMB >= 16 |
|
|| r % 11 == 0 || r % 13 == 0 |
|
#endif |
|
#if BITS_PER_MP_LIMB >= 32 |
|| r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 |
|| r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 |
#if BITS_PER_MP_LIMB == 64 |
#endif |
|
#if BITS_PER_MP_LIMB >= 64 |
|| r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0 |
|| r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0 |
|| r % 47 == 0 || r % 53 == 0 |
|| r % 47 == 0 || r % 53 == 0 |
#endif |
#endif |
Line 81 mpz_probab_prime_p (n, reps) |
|
Line 88 mpz_probab_prime_p (n, reps) |
|
{ |
{ |
return 0; |
return 0; |
} |
} |
|
#endif /* PP */ |
|
|
/* Do more dividing. We collect small primes, using umul_ppmm, until we |
/* 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, |
overflow a single limb. We divide our number by the small primes product, |
Line 95 mpz_probab_prime_p (n, reps) |
|
Line 103 mpz_probab_prime_p (n, reps) |
|
nprimes = 0; |
nprimes = 0; |
p = 1; |
p = 1; |
ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; |
ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2; |
for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2) |
for (q = PP_FIRST_OMITTED; q < ln2; q += 2) |
{ |
{ |
if (isprime (q)) |
if (isprime (q)) |
{ |
{ |
Line 106 mpz_probab_prime_p (n, reps) |
|
Line 114 mpz_probab_prime_p (n, reps) |
|
while (--nprimes >= 0) |
while (--nprimes >= 0) |
if (r % primes[nprimes] == 0) |
if (r % primes[nprimes] == 0) |
{ |
{ |
if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0) |
ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0); |
abort (); |
|
return 0; |
return 0; |
} |
} |
p = q; |
p = q; |
Line 127 mpz_probab_prime_p (n, reps) |
|
Line 134 mpz_probab_prime_p (n, reps) |
|
} |
} |
|
|
static int |
static int |
#if __STDC__ |
|
isprime (unsigned long int t) |
isprime (unsigned long int t) |
#else |
|
isprime (t) |
|
unsigned long int t; |
|
#endif |
|
{ |
{ |
unsigned long int q, r, d; |
unsigned long int q, r, d; |
|
|
|
|
r = t - q * d; |
r = t - q * d; |
if (q < d) |
if (q < d) |
return 1; |
return 1; |
} |
|
return 0; |
|
} |
|
|
|
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__ |
|
mpz_millerrabin (mpz_srcptr n, int reps) |
|
#else |
|
mpz_millerrabin (n, reps) |
|
mpz_srcptr n; |
|
int reps; |
|
#endif |
|
{ |
|
int r; |
|
mpz_t nm1, x, y, q; |
|
unsigned long int k; |
|
gmp_randstate_t rstate; |
|
int is_prime; |
|
TMP_DECL (marker); |
|
TMP_MARK (marker); |
|
|
|
MPZ_TMP_INIT (nm1, SIZ (n) + 1); |
|
mpz_sub_ui (nm1, n, 1L); |
|
|
|
MPZ_TMP_INIT (x, SIZ (n)); |
|
MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */ |
|
|
|
/* Perform a Fermat test. */ |
|
mpz_set_ui (x, 210L); |
|
mpz_powm (y, x, nm1, n); |
|
if (mpz_cmp_ui (y, 1L) != 0) |
|
{ |
|
TMP_FREE (marker); |
|
return 0; |
|
} |
|
|
|
MPZ_TMP_INIT (q, SIZ (n)); |
|
|
|
/* Find q and k, where q is odd and n = 1 + 2**k * q. */ |
|
k = mpz_scan1 (nm1, 0L); |
|
mpz_tdiv_q_2exp (q, nm1, k); |
|
|
|
gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L); |
|
|
|
is_prime = 1; |
|
for (r = 0; r < reps && is_prime; r++) |
|
{ |
|
do |
|
mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1); |
|
while (mpz_cmp_ui (x, 1L) <= 0); |
|
|
|
is_prime = millerrabin (n, nm1, x, y, q, k); |
|
} |
|
|
|
gmp_randclear (rstate); |
|
|
|
TMP_FREE (marker); |
|
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; |
return 0; |
} |
} |