Annotation of OpenXM_contrib/gmp/mpz/pprime_p.c, Revision 1.1.1.4
1.1 maekawa 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'
1.1.1.2 maekawa 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.
1.1 maekawa 8:
1.1.1.4 ! ohara 9: Copyright 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000, 2001 Free Software
1.1.1.2 maekawa 10: Foundation, Inc. Miller-Rabin code contributed by John Amanatides.
1.1 maekawa 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
1.1.1.2 maekawa 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
1.1 maekawa 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
1.1.1.2 maekawa 21: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 22: License for more details.
23:
1.1.1.2 maekawa 24: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 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. */
28:
29: #include "gmp.h"
1.1.1.2 maekawa 30: #include "gmp-impl.h"
31: #include "longlong.h"
1.1 maekawa 32:
1.1.1.2 maekawa 33: static int isprime _PROTO ((unsigned long int t));
34:
35: int
36: mpz_probab_prime_p (mpz_srcptr n, int reps)
1.1 maekawa 37: {
1.1.1.2 maekawa 38: mp_limb_t r;
1.1 maekawa 39:
1.1.1.2 maekawa 40: /* Handle small and negative n. */
41: if (mpz_cmp_ui (n, 1000000L) <= 0)
1.1 maekawa 42: {
1.1.1.2 maekawa 43: int is_prime;
44: if (mpz_sgn (n) < 0)
45: {
46: /* Negative number. Negate and call ourselves. */
47: mpz_t n2;
48: mpz_init (n2);
49: mpz_neg (n2, n);
50: is_prime = mpz_probab_prime_p (n2, reps);
51: mpz_clear (n2);
52: return is_prime;
53: }
54: is_prime = isprime (mpz_get_ui (n));
55: return is_prime ? 2 : 0;
1.1 maekawa 56: }
57:
1.1.1.2 maekawa 58: /* If n is now even, it is not a prime. */
59: if ((mpz_get_ui (n) & 1) == 0)
60: return 0;
1.1 maekawa 61:
1.1.1.4 ! ohara 62: #if defined (PP)
1.1.1.2 maekawa 63: /* Check if n has small factors. */
1.1.1.4 ! ohara 64: #if defined (PP_INVERTED)
! 65: r = MPN_MOD_OR_PREINV_MOD_1 (PTR(n), SIZ(n), (mp_limb_t) PP,
! 66: (mp_limb_t) PP_INVERTED);
! 67: #else
! 68: r = mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) PP);
! 69: #endif
! 70: if (r % 3 == 0
! 71: #if BITS_PER_MP_LIMB >= 4
! 72: || r % 5 == 0
! 73: #endif
! 74: #if BITS_PER_MP_LIMB >= 8
! 75: || r % 7 == 0
! 76: #endif
! 77: #if BITS_PER_MP_LIMB >= 16
! 78: || r % 11 == 0 || r % 13 == 0
! 79: #endif
! 80: #if BITS_PER_MP_LIMB >= 32
1.1.1.2 maekawa 81: || r % 17 == 0 || r % 19 == 0 || r % 23 == 0 || r % 29 == 0
1.1.1.4 ! ohara 82: #endif
! 83: #if BITS_PER_MP_LIMB >= 64
1.1.1.2 maekawa 84: || r % 31 == 0 || r % 37 == 0 || r % 41 == 0 || r % 43 == 0
85: || r % 47 == 0 || r % 53 == 0
86: #endif
87: )
88: {
89: return 0;
90: }
1.1.1.4 ! ohara 91: #endif /* PP */
1.1 maekawa 92:
1.1.1.2 maekawa 93: /* Do more dividing. We collect small primes, using umul_ppmm, until we
94: overflow a single limb. We divide our number by the small primes product,
95: and look for factors in the remainder. */
96: {
97: unsigned long int ln2;
98: unsigned long int q;
99: mp_limb_t p1, p0, p;
100: unsigned int primes[15];
101: int nprimes;
102:
103: nprimes = 0;
104: p = 1;
105: ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
1.1.1.4 ! ohara 106: for (q = PP_FIRST_OMITTED; q < ln2; q += 2)
1.1.1.2 maekawa 107: {
108: if (isprime (q))
109: {
110: umul_ppmm (p1, p0, p, q);
111: if (p1 != 0)
112: {
113: r = mpn_mod_1 (PTR(n), SIZ(n), p);
114: while (--nprimes >= 0)
115: if (r % primes[nprimes] == 0)
116: {
1.1.1.4 ! ohara 117: ASSERT_ALWAYS (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) == 0);
1.1.1.2 maekawa 118: return 0;
119: }
120: p = q;
121: nprimes = 0;
122: }
123: else
124: {
125: p = p0;
126: }
127: primes[nprimes++] = q;
128: }
129: }
130: }
131:
132: /* Perform a number of Miller-Rabin tests. */
133: return mpz_millerrabin (n, reps);
134: }
135:
136: static int
137: isprime (unsigned long int t)
138: {
139: unsigned long int q, r, d;
140:
141: if (t < 3 || (t & 1) == 0)
142: return t == 2;
143:
144: for (d = 3, r = 1; r != 0; d += 2)
1.1 maekawa 145: {
1.1.1.2 maekawa 146: q = t / d;
147: r = t - q * d;
148: if (q < d)
1.1 maekawa 149: return 1;
1.1.1.2 maekawa 150: }
151: return 0;
1.1 maekawa 152: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pprime_p.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz