1 /* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
3 Copyright 1999, 2000, 2001, 2008, 2009 Free Software Foundation, Inc.
5 Contributed to the GNU project by Niels Möller and Torbjörn Granlund.
7 This file is part of the GNU MP Library.
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MP Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
26 static const unsigned char primegap[] =
28 2,2,4,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6,
29 2,10,2,6,6,4,6,6,2,10,2,4,2,12,12,4,2,4,6,2,10,6,6,6,2,6,4,2,10,14,4,2,
30 4,14,6,10,2,4,6,8,6,6,4,6,8,4,8,10,2,10,2,6,4,6,8,4,2,4,12,8,4,8,4,6,
31 12,2,18,6,10,6,6,2,6,10,6,6,2,6,6,4,2,12,10,2,4,6,6,2,12,4,6,8,10,8,10,8,
32 6,6,4,8,6,4,8,4,14,10,12,2,10,2,4,2,10,14,4,2,4,14,4,2,4,20,4,8,10,8,4,6,
36 #define NUMBER_OF_PRIMES 167
39 mpz_nextprime (mpz_ptr p, mpz_srcptr n)
41 unsigned short *moduli;
42 unsigned long difference;
52 /* First handle tiny numbers */
53 if (mpz_cmp_ui (n, 2) < 0)
61 if (mpz_cmp_ui (p, 7) <= 0)
65 count_leading_zeros (cnt, PTR(p)[pn - 1]);
66 nbits = pn * GMP_NUMB_BITS - (cnt - GMP_NAIL_BITS);
67 if (nbits / 2 >= NUMBER_OF_PRIMES)
68 prime_limit = NUMBER_OF_PRIMES - 1;
70 prime_limit = nbits / 2;
74 /* Compute residues modulo small odd primes */
75 moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short);
79 /* FIXME: Compute lazily? */
81 for (i = 0; i < prime_limit; i++)
83 moduli[i] = mpz_fdiv_ui (p, prime);
87 #define INCR_LIMIT 0x10000 /* deep science */
89 for (difference = incr = 0; incr < INCR_LIMIT; difference += 2)
91 /* First check residues */
93 for (i = 0; i < prime_limit; i++)
96 /* FIXME: Reduce moduli + incr and store back, to allow for
97 division-free reductions. Alternatively, table primes[]'s
98 inverses (mod 2^16). */
99 r = (moduli[i] + incr) % prime;
100 prime += primegap[i];
106 mpz_add_ui (p, p, difference);
109 /* Miller-Rabin test */
110 if (mpz_millerrabin (p, 10))
115 mpz_add_ui (p, p, difference);