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 GMP_LIMB_BITS >= 4
76 #if GMP_LIMB_BITS >= 8
79 #if GMP_LIMB_BITS >= 16
80 || r
% 11 == 0 || r
% 13 == 0
82 #if GMP_LIMB_BITS >= 32
83 || r
% 17 == 0 || r
% 19 == 0 || r
% 23 == 0 || r
% 29 == 0
85 #if GMP_LIMB_BITS >= 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)