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. 8 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. 11 12 This file is part of the GNU MP Library. 13 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. 18 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. 23 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/. */ 26 27 #include "gmp.h" 28 #include "gmp-impl.h" 29 #include "longlong.h" 30 31 static int isprime __GMP_PROTO ((unsigned long int)); 32 33 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. */ 38 39 int 40 mpz_probab_prime_p (mpz_srcptr n, int reps) 41 { 42 mp_limb_t r; 43 mpz_t n2; 44 45 /* Handle small and negative n. */ 46 if (mpz_cmp_ui (n, 1000000L) <= 0) 47 { 48 int is_prime; 49 if (mpz_cmpabs_ui (n, 1000000L) <= 0) 50 { 51 is_prime = isprime (mpz_get_ui (n)); 52 return is_prime ? 2 : 0; 53 } 54 /* Negative number. Negate and fall out. */ 55 PTR(n2) = PTR(n); 56 SIZ(n2) = -SIZ(n); 57 n = n2; 58 } 59 60 /* If n is now even, it is not a prime. */ 61 if ((mpz_get_ui (n) & 1) == 0) 62 return 0; 63 64 #if defined (PP) 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); 69 #else 70 r = mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) PP); 71 #endif 72 if (r % 3 == 0 73 #if GMP_LIMB_BITS >= 4 74 || r % 5 == 0 75 #endif 76 #if GMP_LIMB_BITS >= 8 77 || r % 7 == 0 78 #endif 79 #if GMP_LIMB_BITS >= 16 80 || r % 11 == 0 || r % 13 == 0 81 #endif 82 #if GMP_LIMB_BITS >= 32 83 || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0 84 #endif 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 88 #endif 89 ) 90 { 91 return 0; 92 } 93 #endif /* PP */ 94 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. */ 98 { 99 unsigned long int ln2; 100 unsigned long int q; 101 mp_limb_t p1, p0, p; 102 unsigned int primes[15]; 103 int nprimes; 104 105 nprimes = 0; 106 p = 1; 107 ln2 = mpz_sizeinbase (n, 2); /* FIXME: tune this limit */ 108 for (q = PP_FIRST_OMITTED; q < ln2; q += 2) 109 { 110 if (isprime (q)) 111 { 112 umul_ppmm (p1, p0, p, q); 113 if (p1 != 0) 114 { 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) 118 { 119 ASSERT_ALWAYS (mpn_mod_1 (PTR(n), (mp_size_t) SIZ(n), (mp_limb_t) primes[nprimes]) == 0); 120 return 0; 121 } 122 p = q; 123 nprimes = 0; 124 } 125 else 126 { 127 p = p0; 128 } 129 primes[nprimes++] = q; 130 } 131 } 132 } 133 134 /* Perform a number of Miller-Rabin tests. */ 135 return mpz_millerrabin (n, reps); 136 } 137 138 static int 139 isprime (unsigned long int t) 140 { 141 unsigned long int q, r, d; 142 143 if (t < 3 || (t & 1) == 0) 144 return t == 2; 145 146 for (d = 3, r = 1; r != 0; d += 2) 147 { 148 q = t / d; 149 r = t - q * d; 150 if (q < d) 151 return 1; 152 } 153 return 0; 154 } 155