Annotation of OpenXM/src/kan96xx/gmp-2.0.2/mpz/pprime_p.c, Revision 1.1.1.1
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'
5: prime. The probability of a false positive is (1/4)**reps, where
6: reps is the number of internal passes of the probabilistic algorithm.
7: Knuth indicates that 25 passes are reasonable.
8:
9: Copyright (C) 1991, 1993, 1994 Free Software Foundation, Inc.
10: 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 Library General Public License as published by
16: the Free Software Foundation; either version 2 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 Library General Public
22: License for more details.
23:
24: You should have received a copy of the GNU Library 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. */
28:
29: #include "gmp.h"
30:
31: static int
32: possibly_prime (n, n_minus_1, x, y, q, k)
33: mpz_srcptr n;
34: mpz_srcptr n_minus_1;
35: mpz_ptr x;
36: mpz_ptr y;
37: mpz_srcptr q;
38: unsigned long int k;
39: {
40: unsigned long int i;
41:
42: /* find random x s.t. 1 < x < n */
43: do
44: {
45: mpz_random (x, mpz_size (n));
46: mpz_mmod (x, x, n);
47: }
48: while (mpz_cmp_ui (x, 1L) <= 0);
49:
50: mpz_powm (y, x, q, n);
51:
52: if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, n_minus_1) == 0)
53: return 1;
54:
55: for (i = 1; i < k; i++)
56: {
57: mpz_powm_ui (y, y, 2L, n);
58: if (mpz_cmp (y, n_minus_1) == 0)
59: return 1;
60: if (mpz_cmp_ui (y, 1L) == 0)
61: return 0;
62: }
63: return 0;
64: }
65:
66: int
67: #if __STDC__
68: mpz_probab_prime_p (mpz_srcptr m, int reps)
69: #else
70: mpz_probab_prime_p (m, reps)
71: mpz_srcptr m;
72: int reps;
73: #endif
74: {
75: mpz_t n, n_minus_1, x, y, q;
76: int i, is_prime;
77: unsigned long int k;
78:
79: mpz_init (n);
80: /* Take the absolute value of M, to handle positive and negative primes. */
81: mpz_abs (n, m);
82:
83: if (mpz_cmp_ui (n, 3L) <= 0)
84: {
85: mpz_clear (n);
86: return mpz_cmp_ui (n, 1L) > 0;
87: }
88:
89: if ((mpz_get_ui (n) & 1) == 0)
90: {
91: mpz_clear (n);
92: return 0; /* even */
93: }
94:
95: mpz_init (n_minus_1);
96: mpz_sub_ui (n_minus_1, n, 1L);
97: mpz_init (x);
98: mpz_init (y);
99:
100: /* find q and k, s.t. n = 1 + 2**k * q */
101: mpz_init_set (q, n_minus_1);
102: k = mpz_scan1 (q, 0);
103: mpz_tdiv_q_2exp (q, q, k);
104:
105: is_prime = 1;
106: for (i = 0; i < reps && is_prime; i++)
107: is_prime &= possibly_prime (n, n_minus_1, x, y, q, k);
108:
109: mpz_clear (n_minus_1);
110: mpz_clear (n);
111: mpz_clear (x);
112: mpz_clear (y);
113: mpz_clear (q);
114: return is_prime;
115: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>