1/* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
2
3Copyright 1999, 2000, 2001, 2008, 2009 Free Software Foundation, Inc.
4
5Contributed to the GNU project by Niels M�ller and Torbjorn Granlund.
6
7This file is part of the GNU MP Library.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MP Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
21
22#include "gmp.h"
23#include "gmp-impl.h"
24#include "longlong.h"
25
26static const unsigned char primegap[] =
27{
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,
33  6,14,4,6,6,8,6,12
34};
35
36#define NUMBER_OF_PRIMES 167
37
38void
39mpz_nextprime (mpz_ptr p, mpz_srcptr n)
40{
41  unsigned short *moduli;
42  unsigned long difference;
43  int i;
44  unsigned prime_limit;
45  unsigned long prime;
46  int cnt;
47  mp_size_t pn;
48  mp_bitcnt_t nbits;
49  unsigned incr;
50  TMP_SDECL;
51
52  /* First handle tiny numbers */
53  if (mpz_cmp_ui (n, 2) < 0)
54    {
55      mpz_set_ui (p, 2);
56      return;
57    }
58  mpz_add_ui (p, n, 1);
59  mpz_setbit (p, 0);
60
61  if (mpz_cmp_ui (p, 7) <= 0)
62    return;
63
64  pn = SIZ(p);
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;
69  else
70    prime_limit = nbits / 2;
71
72  TMP_SMARK;
73
74  /* Compute residues modulo small odd primes */
75  moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short);
76
77  for (;;)
78    {
79      /* FIXME: Compute lazily? */
80      prime = 3;
81      for (i = 0; i < prime_limit; i++)
82	{
83	  moduli[i] = mpz_fdiv_ui (p, prime);
84	  prime += primegap[i];
85	}
86
87#define INCR_LIMIT 0x10000	/* deep science */
88
89      for (difference = incr = 0; incr < INCR_LIMIT; difference += 2)
90	{
91	  /* First check residues */
92	  prime = 3;
93	  for (i = 0; i < prime_limit; i++)
94	    {
95	      unsigned r;
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];
101
102	      if (r == 0)
103		goto next;
104	    }
105
106	  mpz_add_ui (p, p, difference);
107	  difference = 0;
108
109	  /* Miller-Rabin test */
110	  if (mpz_millerrabin (p, 25))
111	    goto done;
112	next:;
113	  incr += 2;
114	}
115      mpz_add_ui (p, p, difference);
116      difference = 0;
117    }
118 done:
119  TMP_SFREE;
120}
121