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.
9 Copyright (C) 1991 Free Software Foundation, Inc.
10 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 General Public License as published by
16 the Free Software Foundation; either version 2, or (at your option)
19 The GNU MP Library is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU General Public License for more details.
24 You should have received a copy of the GNU General Public License
25 along with the GNU MP Library; see the file COPYING. If not, write to
26 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
32 possibly_prime (n, n_minus_1, x, y, q, k)
33 MP_INT *n, *n_minus_1, *x, *y, *q;
38 /* find random x s.t. 1 < x < n */
41 mpz_random (x, mpz_size (n));
44 while (mpz_cmp_ui (x, 1) <= 0);
46 mpz_powm (y, x, q, n);
48 if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, n_minus_1) == 0)
51 for (i = 1; i < k; i++)
53 mpz_powm_ui (y, y, 2, n);
54 if (mpz_cmp (y, n_minus_1) == 0)
56 if (mpz_cmp_ui (y, 1) == 0)
63 mpz_probab_prime_p (m, reps)
67 MP_INT n, n_minus_1, x, y, q;
71 /* Take the absolute value of M, to handle positive and negative primes. */
74 if (mpz_cmp_ui (&n, 3) <= 0)
76 if (mpz_cmp_ui (&n, 1) <= 0)
77 return 0; /* smallest prime is 2 */
81 if ((mpz_get_ui (&n) & 1) == 0)
84 mpz_init (&n_minus_1);
85 mpz_sub_ui (&n_minus_1, &n, 1);
89 /* find q and k, s.t. n = 1 + 2**k * q */
90 mpz_init_set (&q, &n_minus_1);
92 while ((mpz_get_ui (&q) & 1) == 0)
95 mpz_div_2exp (&q, &q, 1);
99 for (i = 0; i < reps && is_prime; i++)
100 is_prime &= possibly_prime (&n, &n_minus_1, &x, &y, &q, k);
102 mpz_clear (&n_minus_1);