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. If it returns 2, n is surely prime. The probability of a false
6 positive is (1/4)**reps, where reps is the number of internal passes of the
7 probabilistic algorithm. Knuth indicates that 25 passes are reasonable.
9 Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2005 Free
10 Software Foundation, Inc. Miller-Rabin code contributed by John Amanatides.
12 This file is part of the GNU MP Library.
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
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 Lesser General Public
22 License for more details.
24 You should have received a copy of the GNU Lesser General Public License
25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
31 static int isprime __GMP_PROTO ((unsigned long int));
34 /* MPN_MOD_OR_MODEXACT_1_ODD can be used instead of mpn_mod_1 for the trial
35 division. It gives a result which is not the actual remainder r but a
36 value congruent to r*2^n mod d. Since all the primes being tested are
37 odd, r*2^n mod p will be 0 if and only if r mod p is 0. */
40 mpz_probab_prime_p (mpz_srcptr n, int reps)
45 /* Handle small and negative n. */
46 if (mpz_cmp_ui (n, 1000000L) <= 0)
49 if (mpz_cmpabs_ui (n, 1000000L) <= 0)
51 is_prime = isprime (mpz_get_ui (n));
52 return is_prime ? 2 : 0;
54 /* Negative number. Negate and fall out. */
60 /* If n is now even, it is not a prime. */
61 if ((mpz_get_ui (n) & 1) == 0)
65 /* Check if n has small factors. */
66 #if defined (PP_INVERTED)
67 r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP,
68 (mp_limb_t) PP_INVERTED);
70 r = mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP);
73 #if BITS_PER_MP_LIMB >= 4
76 #if BITS_PER_MP_LIMB >= 8
79 #if BITS_PER_MP_LIMB >= 16
80 || r % 11 == 0 || r % 13 == 0
82 #if BITS_PER_MP_LIMB >= 32
83 || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
85 #if BITS_PER_MP_LIMB >= 64
86 || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
87 || r % 47 == 0 || r % 53 == 0
95 /* Do more dividing. We collect small primes, using umul_ppmm, until we
96 overflow a single limb. We divide our number by the small primes product,
97 and look for factors in the remainder. */
99 unsigned long int ln2;
102 unsigned int primes[15];
107 ln2 = mpz_sizeinbase (n, 2); /* FIXME: tune this limit */
108 for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
112 umul_ppmm (p1, p0, p, q);
115 r = MPN_MOD_OR_MODEXACT_1_ODD (PTR(n), (mp_size_t) SIZ(n), p);
116 while (--nprimes >= 0)
117 if (r % primes[nprimes] == 0)
119 ASSERT_ALWAYS (mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
129 primes[nprimes++] = q;
134 /* Perform a number of Miller-Rabin tests. */
135 return mpz_millerrabin (n, reps);
139 isprime (unsigned long int t)
141 unsigned long int q, r, d;
143 if (t < 3 || (t & 1) == 0)
146 for (d = 3, r = 1; r != 0; d += 2)