Annotation of OpenXM_contrib/gmp/mpz/pprime_p.c, Revision 1.1.1.2
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: return 0;
! 185: TMP_FREE (marker);
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