|
Packit |
5c3484 |
/* mpz_probab_prime_p --
|
|
Packit |
5c3484 |
An implementation of the probabilistic primality test found in Knuth's
|
|
Packit |
5c3484 |
Seminumerical Algorithms book. If the function mpz_probab_prime_p()
|
|
Packit |
5c3484 |
returns 0 then n is not prime. If it returns 1, then n is 'probably'
|
|
Packit |
5c3484 |
prime. If it returns 2, n is surely prime. The probability of a false
|
|
Packit |
5c3484 |
positive is (1/4)**reps, where reps is the number of internal passes of the
|
|
Packit |
5c3484 |
probabilistic algorithm. Knuth indicates that 25 passes are reasonable.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
Copyright 1991, 1993, 1994, 1996-2002, 2005 Free Software Foundation, Inc.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
This file is part of the GNU MP Library.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is free software; you can redistribute it and/or modify
|
|
Packit |
5c3484 |
it under the terms of either:
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU Lesser General Public License as published by the Free
|
|
Packit |
5c3484 |
Software Foundation; either version 3 of the License, or (at your
|
|
Packit |
5c3484 |
option) any later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
* the GNU General Public License as published by the Free Software
|
|
Packit |
5c3484 |
Foundation; either version 2 of the License, or (at your option) any
|
|
Packit |
5c3484 |
later version.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
or both in parallel, as here.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
The GNU MP Library is distributed in the hope that it will be useful, but
|
|
Packit |
5c3484 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
Packit |
5c3484 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
5c3484 |
for more details.
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
You should have received copies of the GNU General Public License and the
|
|
Packit |
5c3484 |
GNU Lesser General Public License along with the GNU MP Library. If not,
|
|
Packit |
5c3484 |
see https://www.gnu.org/licenses/. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#include "gmp.h"
|
|
Packit |
5c3484 |
#include "gmp-impl.h"
|
|
Packit |
5c3484 |
#include "longlong.h"
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
static int isprime (unsigned long int);
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* MPN_MOD_OR_MODEXACT_1_ODD can be used instead of mpn_mod_1 for the trial
|
|
Packit |
5c3484 |
division. It gives a result which is not the actual remainder r but a
|
|
Packit |
5c3484 |
value congruent to r*2^n mod d. Since all the primes being tested are
|
|
Packit |
5c3484 |
odd, r*2^n mod p will be 0 if and only if r mod p is 0. */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
int
|
|
Packit |
5c3484 |
mpz_probab_prime_p (mpz_srcptr n, int reps)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
mp_limb_t r;
|
|
Packit |
5c3484 |
mpz_t n2;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Handle small and negative n. */
|
|
Packit |
5c3484 |
if (mpz_cmp_ui (n, 1000000L) <= 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
int is_prime;
|
|
Packit |
5c3484 |
if (mpz_cmpabs_ui (n, 1000000L) <= 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
is_prime = isprime (mpz_get_ui (n));
|
|
Packit |
5c3484 |
return is_prime ? 2 : 0;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
/* Negative number. Negate and fall out. */
|
|
Packit |
5c3484 |
PTR(n2) = PTR(n);
|
|
Packit |
5c3484 |
SIZ(n2) = -SIZ(n);
|
|
Packit |
5c3484 |
n = n2;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* If n is now even, it is not a prime. */
|
|
Packit |
5c3484 |
if ((mpz_get_ui (n) & 1) == 0)
|
|
Packit |
5c3484 |
return 0;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
#if defined (PP)
|
|
Packit |
5c3484 |
/* Check if n has small factors. */
|
|
Packit |
5c3484 |
#if defined (PP_INVERTED)
|
|
Packit |
5c3484 |
r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP,
|
|
Packit |
5c3484 |
(mp_limb_t) PP_INVERTED);
|
|
Packit |
5c3484 |
#else
|
|
Packit |
5c3484 |
r = mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP);
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
if (r % 3 == 0
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS >= 4
|
|
Packit |
5c3484 |
|| r % 5 == 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS >= 8
|
|
Packit |
5c3484 |
|| r % 7 == 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS >= 16
|
|
Packit |
5c3484 |
|| r % 11 == 0 || r % 13 == 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS >= 32
|
|
Packit |
5c3484 |
|| r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
#if GMP_LIMB_BITS >= 64
|
|
Packit |
5c3484 |
|| r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
|
|
Packit |
5c3484 |
|| r % 47 == 0 || r % 53 == 0
|
|
Packit |
5c3484 |
#endif
|
|
Packit |
5c3484 |
)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
return 0;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
#endif /* PP */
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Do more dividing. We collect small primes, using umul_ppmm, until we
|
|
Packit |
5c3484 |
overflow a single limb. We divide our number by the small primes product,
|
|
Packit |
5c3484 |
and look for factors in the remainder. */
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
unsigned long int ln2;
|
|
Packit |
5c3484 |
unsigned long int q;
|
|
Packit |
5c3484 |
mp_limb_t p1, p0, p;
|
|
Packit |
5c3484 |
unsigned int primes[15];
|
|
Packit |
5c3484 |
int nprimes;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
nprimes = 0;
|
|
Packit |
5c3484 |
p = 1;
|
|
Packit |
5c3484 |
ln2 = mpz_sizeinbase (n, 2); /* FIXME: tune this limit */
|
|
Packit |
5c3484 |
for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
if (isprime (q))
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
umul_ppmm (p1, p0, p, q);
|
|
Packit |
5c3484 |
if (p1 != 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
r = MPN_MOD_OR_MODEXACT_1_ODD (PTR(n), (mp_size_t) SIZ(n), p);
|
|
Packit |
5c3484 |
while (--nprimes >= 0)
|
|
Packit |
5c3484 |
if (r % primes[nprimes] == 0)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
ASSERT_ALWAYS (mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
|
|
Packit |
5c3484 |
return 0;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
p = q;
|
|
Packit |
5c3484 |
nprimes = 0;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
else
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
p = p0;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
primes[nprimes++] = q;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
/* Perform a number of Miller-Rabin tests. */
|
|
Packit |
5c3484 |
return mpz_millerrabin (n, reps);
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
static int
|
|
Packit |
5c3484 |
isprime (unsigned long int t)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
unsigned long int q, r, d;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
if (t < 3 || (t & 1) == 0)
|
|
Packit |
5c3484 |
return t == 2;
|
|
Packit |
5c3484 |
|
|
Packit |
5c3484 |
for (d = 3, r = 1; r != 0; d += 2)
|
|
Packit |
5c3484 |
{
|
|
Packit |
5c3484 |
q = t / d;
|
|
Packit |
5c3484 |
r = t - q * d;
|
|
Packit |
5c3484 |
if (q < d)
|
|
Packit |
5c3484 |
return 1;
|
|
Packit |
5c3484 |
}
|
|
Packit |
5c3484 |
return 0;
|
|
Packit |
5c3484 |
}
|