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 (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software
10 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 2.1 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; see the file COPYING.LIB. If not, write to
26 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27 MA 02111-1307, USA. */
33 static int isprime _PROTO ((unsigned long int t));
34 static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps));
38 mpz_probab_prime_p (mpz_srcptr n, int reps)
40 mpz_probab_prime_p (n, reps)
47 /* Handle small and negative n. */
48 if (mpz_cmp_ui (n, 1000000L) <= 0)
53 /* Negative number. Negate and call ourselves. */
57 is_prime = mpz_probab_prime_p (n2, reps);
61 is_prime = isprime (mpz_get_ui (n));
62 return is_prime ? 2 : 0;
65 /* If n is now even, it is not a prime. */
66 if ((mpz_get_ui (n) & 1) == 0)
69 /* Check if n has small factors. */
70 if (UDIV_TIME > (2 * UMUL_TIME + 6))
71 r = mpn_preinv_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP, (mp_limb_t) PP_INVERTED);
73 r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
74 if (r % 3 == 0 || r % 5 == 0 || r % 7 == 0 || r % 11 == 0 || r % 13 == 0
75 || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
76 #if BITS_PER_MP_LIMB == 64
77 || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
78 || r % 47 == 0 || r % 53 == 0
85 /* Do more dividing. We collect small primes, using umul_ppmm, until we
86 overflow a single limb. We divide our number by the small primes product,
87 and look for factors in the remainder. */
89 unsigned long int ln2;
92 unsigned int primes[15];
97 ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
98 for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2)
102 umul_ppmm (p1, p0, p, q);
105 r = mpn_mod_1 (PTR(n), SIZ(n), p);
106 while (--nprimes >= 0)
107 if (r % primes[nprimes] == 0)
109 if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0)
120 primes[nprimes++] = q;
125 /* Perform a number of Miller-Rabin tests. */
126 return mpz_millerrabin (n, reps);
131 isprime (unsigned long int t)
137 unsigned long int q, r, d;
139 if (t < 3 || (t & 1) == 0)
142 for (d = 3, r = 1; r != 0; d += 2)
152 static int millerrabin _PROTO ((mpz_srcptr n, mpz_srcptr nm1,
153 mpz_ptr x, mpz_ptr y,
154 mpz_srcptr q, unsigned long int k));
158 mpz_millerrabin (mpz_srcptr n, int reps)
160 mpz_millerrabin (n, reps)
168 gmp_randstate_t rstate;
173 MPZ_TMP_INIT (nm1, SIZ (n) + 1);
174 mpz_sub_ui (nm1, n, 1L);
176 MPZ_TMP_INIT (x, SIZ (n));
177 MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
179 /* Perform a Fermat test. */
180 mpz_set_ui (x, 210L);
181 mpz_powm (y, x, nm1, n);
182 if (mpz_cmp_ui (y, 1L) != 0)
188 MPZ_TMP_INIT (q, SIZ (n));
190 /* Find q and k, where q is odd and n = 1 + 2**k * q. */
191 k = mpz_scan1 (nm1, 0L);
192 mpz_tdiv_q_2exp (q, nm1, k);
194 gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L);
197 for (r = 0; r < reps && is_prime; r++)
200 mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1);
201 while (mpz_cmp_ui (x, 1L) <= 0);
203 is_prime = millerrabin (n, nm1, x, y, q, k);
206 gmp_randclear (rstate);
214 millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
215 mpz_srcptr q, unsigned long int k)
217 millerrabin (n, nm1, x, y, q, k)
228 mpz_powm (y, x, q, n);
230 if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
233 for (i = 1; i < k; i++)
235 mpz_powm_ui (y, y, 2L, n);
236 if (mpz_cmp (y, nm1) == 0)
238 if (mpz_cmp_ui (y, 1L) == 0)