Annotation of OpenXM_contrib/gmp/mpz/pprime_p.c, Revision 1.1.1.3
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.2 maekawa 9: Copyright (C) 1991, 1993, 1994, 1996, 1997, 1998, 1999, 2000 Free Software
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: static int mpz_millerrabin _PROTO ((mpz_srcptr n, int reps));
35:
36: int
37: #if __STDC__
38: mpz_probab_prime_p (mpz_srcptr n, int reps)
39: #else
40: mpz_probab_prime_p (n, reps)
1.1 maekawa 41: mpz_srcptr n;
1.1.1.2 maekawa 42: int reps;
43: #endif
1.1 maekawa 44: {
1.1.1.2 maekawa 45: mp_limb_t r;
1.1 maekawa 46:
1.1.1.2 maekawa 47: /* Handle small and negative n. */
48: if (mpz_cmp_ui (n, 1000000L) <= 0)
1.1 maekawa 49: {
1.1.1.2 maekawa 50: int is_prime;
51: if (mpz_sgn (n) < 0)
52: {
53: /* Negative number. Negate and call ourselves. */
54: mpz_t n2;
55: mpz_init (n2);
56: mpz_neg (n2, n);
57: is_prime = mpz_probab_prime_p (n2, reps);
58: mpz_clear (n2);
59: return is_prime;
60: }
61: is_prime = isprime (mpz_get_ui (n));
62: return is_prime ? 2 : 0;
1.1 maekawa 63: }
64:
1.1.1.2 maekawa 65: /* If n is now even, it is not a prime. */
66: if ((mpz_get_ui (n) & 1) == 0)
67: return 0;
1.1 maekawa 68:
1.1.1.2 maekawa 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);
72: else
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
79: #endif
80: )
81: {
82: return 0;
83: }
1.1 maekawa 84:
1.1.1.2 maekawa 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. */
88: {
89: unsigned long int ln2;
90: unsigned long int q;
91: mp_limb_t p1, p0, p;
92: unsigned int primes[15];
93: int nprimes;
94:
95: nprimes = 0;
96: p = 1;
97: ln2 = mpz_sizeinbase (n, 2) / 30; ln2 = ln2 * ln2;
98: for (q = BITS_PER_MP_LIMB == 64 ? 59 : 31; q < ln2; q += 2)
99: {
100: if (isprime (q))
101: {
102: umul_ppmm (p1, p0, p, q);
103: if (p1 != 0)
104: {
105: r = mpn_mod_1 (PTR(n), SIZ(n), p);
106: while (--nprimes >= 0)
107: if (r % primes[nprimes] == 0)
108: {
109: if (mpn_mod_1 (PTR(n), SIZ(n), (mp_limb_t) primes[nprimes]) != 0)
110: abort ();
111: return 0;
112: }
113: p = q;
114: nprimes = 0;
115: }
116: else
117: {
118: p = p0;
119: }
120: primes[nprimes++] = q;
121: }
122: }
123: }
124:
125: /* Perform a number of Miller-Rabin tests. */
126: return mpz_millerrabin (n, reps);
127: }
128:
129: static int
130: #if __STDC__
131: isprime (unsigned long int t)
132: #else
133: isprime (t)
134: unsigned long int t;
135: #endif
136: {
137: unsigned long int q, r, d;
138:
139: if (t < 3 || (t & 1) == 0)
140: return t == 2;
141:
142: for (d = 3, r = 1; r != 0; d += 2)
1.1 maekawa 143: {
1.1.1.2 maekawa 144: q = t / d;
145: r = t - q * d;
146: if (q < d)
1.1 maekawa 147: return 1;
148: }
149: return 0;
150: }
151:
1.1.1.2 maekawa 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));
155:
156: static int
1.1 maekawa 157: #if __STDC__
1.1.1.2 maekawa 158: mpz_millerrabin (mpz_srcptr n, int reps)
1.1 maekawa 159: #else
1.1.1.2 maekawa 160: mpz_millerrabin (n, reps)
161: mpz_srcptr n;
1.1 maekawa 162: int reps;
163: #endif
164: {
1.1.1.2 maekawa 165: int r;
166: mpz_t nm1, x, y, q;
1.1 maekawa 167: unsigned long int k;
1.1.1.2 maekawa 168: gmp_randstate_t rstate;
169: int is_prime;
170: TMP_DECL (marker);
171: TMP_MARK (marker);
172:
173: MPZ_TMP_INIT (nm1, SIZ (n) + 1);
174: mpz_sub_ui (nm1, n, 1L);
175:
176: MPZ_TMP_INIT (x, SIZ (n));
177: MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
178:
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)
1.1 maekawa 183: {
1.1.1.2 maekawa 184: TMP_FREE (marker);
1.1.1.3 ! maekawa 185: return 0;
1.1 maekawa 186: }
187:
1.1.1.2 maekawa 188: MPZ_TMP_INIT (q, SIZ (n));
1.1 maekawa 189:
1.1.1.2 maekawa 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);
193:
194: gmp_randinit (rstate, GMP_RAND_ALG_DEFAULT, 32L);
1.1 maekawa 195:
196: is_prime = 1;
1.1.1.2 maekawa 197: for (r = 0; r < reps && is_prime; r++)
198: {
199: do
200: mpz_urandomb (x, rstate, mpz_sizeinbase (n, 2) - 1);
201: while (mpz_cmp_ui (x, 1L) <= 0);
202:
203: is_prime = millerrabin (n, nm1, x, y, q, k);
204: }
1.1 maekawa 205:
1.1.1.2 maekawa 206: gmp_randclear (rstate);
207:
208: TMP_FREE (marker);
1.1 maekawa 209: return is_prime;
1.1.1.2 maekawa 210: }
211:
212: static int
213: #if __STDC__
214: millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
215: mpz_srcptr q, unsigned long int k)
216: #else
217: millerrabin (n, nm1, x, y, q, k)
218: mpz_srcptr n;
219: mpz_srcptr nm1;
220: mpz_ptr x;
221: mpz_ptr y;
222: mpz_srcptr q;
223: unsigned long int k;
224: #endif
225: {
226: unsigned long int i;
227:
228: mpz_powm (y, x, q, n);
229:
230: if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
231: return 1;
232:
233: for (i = 1; i < k; i++)
234: {
235: mpz_powm_ui (y, y, 2L, n);
236: if (mpz_cmp (y, nm1) == 0)
237: return 1;
238: if (mpz_cmp_ui (y, 1L) == 0)
239: return 0;
240: }
241: return 0;
1.1 maekawa 242: }
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