1 /* mpz_nextprime(p,t) - compute the next prime > t and store that in p.
2
3 Copyright 1999-2001, 2008, 2009, 2012, 2020-2022 Free Software
4 Foundation, Inc.
5
6 Contributed to the GNU project by Niels Möller and Torbjorn Granlund.
7 Improved by Seth Troisi.
8
9 This file is part of the GNU MP Library.
10
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of either:
13
14 * the GNU Lesser General Public License as published by the Free
15 Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 or
19
20 * the GNU General Public License as published by the Free Software
21 Foundation; either version 2 of the License, or (at your option) any
22 later version.
23
24 or both in parallel, as here.
25
26 The GNU MP Library is distributed in the hope that it will be useful, but
27 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29 for more details.
30
31 You should have received copies of the GNU General Public License and the
32 GNU Lesser General Public License along with the GNU MP Library. If not,
33 see https://www.gnu.org/licenses/. */
34
35 #include <string.h>
36
37 #include "gmp-impl.h"
38 #include "longlong.h"
39
40 /*********************************************************/
41 /* Section sieve: sieving functions and tools for primes */
42 /*********************************************************/
43
44 static mp_limb_t
45 n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
46
47 static mp_size_t
48 primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
49
50
51 static const unsigned char primegap_small[] =
52 {
53 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,
54 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,
55 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,
56 12,2,18,6,10
57 };
58
59 #define NUMBER_OF_PRIMES 100
60 #define LAST_PRIME 557
61 /* NP_SMALL_LIMIT = prevprime (LAST_PRIME ^ 2) */
62 #define NP_SMALL_LIMIT 310243
63
64 static unsigned long
65 calculate_sievelimit(mp_bitcnt_t nbits) {
66 unsigned long sieve_limit;
67
68 /* Estimate a good sieve bound. Based on derivative of
69 * Merten's 3rd theorem * avg gap * cost of mod
70 * vs
71 * Cost of PRP test O(N^2.55)
72 */
73 if (nbits < 12818)
74 {
75 mpz_t tmp;
76 /* sieve_limit ~= nbits ^ (5/2) / 124 */
77 mpz_init (tmp);
78 mpz_ui_pow_ui (tmp, nbits, 5);
79 mpz_tdiv_q_ui(tmp, tmp, 124*124);
80 /* tmp < 12818^5/(124*124) < 2^55 < 2^64 */
81 mpz_sqrt (tmp, tmp);
82
83 sieve_limit = mpz_get_ui(tmp);
84 mpz_clear (tmp);
85 }
86 else
87 {
88 /* Larger threshold is faster but takes (n/ln(n) + n/24) memory.
89 * For 33,000 bits limitting to 150M is ~12% slower than using the
90 * optimal 1.5G sieve_limit.
91 */
92 sieve_limit = 150000001;
93 }
94
95 ASSERT (1000 < sieve_limit && sieve_limit <= 150000001);
96 return sieve_limit;
97 }
98
99 static unsigned
100 findnext_small (unsigned t, short diff)
101 {
102 /* For diff= 2, expect t = 1 if operand was negative.
103 * For diff=-2, expect t >= 3
104 */
105 ASSERT (t >= 3 || (diff > 0 && t >= 1));
106 ASSERT (t < NP_SMALL_LIMIT);
107
108 /* Start from next candidate (2 or odd) */
109 t = diff > 0 ?
110 (t + 1) | (t != 1) :
111 ((t - 2) | 1) + (t == 3);
112
113 for (; ; t += diff)
114 {
115 unsigned prime = 3;
116 for (int i = 0; ; prime += primegap_small[i++])
117 {
118 unsigned q, r;
119 q = t / prime;
120 r = t - q * prime; /* r = t % prime; */
121 if (q < prime)
122 return t;
123 if (r == 0)
124 break;
125 ASSERT (i < NUMBER_OF_PRIMES);
126 }
127 }
128 }
129
130 static int
131 findnext (mpz_ptr p,
132 unsigned long(*negative_mod_ui)(const mpz_t, unsigned long),
133 void(*increment_ui)(mpz_t, const mpz_t, unsigned long))
134 {
135 char *composite;
136 const unsigned char *primegap;
137 unsigned long prime_limit;
138 mp_size_t pn;
139 mp_bitcnt_t nbits;
140 int i, m;
141 unsigned odds_in_composite_sieve;
142 TMP_DECL;
143
144 TMP_MARK;
145 pn = SIZ(p);
146 MPN_SIZEINBASE_2EXP(nbits, PTR(p), pn, 1);
147 /* Smaller numbers handled earlier */
148 ASSERT (nbits >= 3);
149 /* Make p odd */
150 PTR(p)[0] |= 1;
151
152 if (nbits / 2 <= NUMBER_OF_PRIMES)
153 {
154 primegap = primegap_small;
155 prime_limit = nbits / 2;
156 }
157 else
158 {
159 unsigned long sieve_limit;
160 mp_limb_t *sieve;
161 unsigned char *primegap_tmp;
162 unsigned long last_prime;
163
164 /* sieve numbers up to sieve_limit and save prime count */
165 sieve_limit = calculate_sievelimit(nbits);
166 sieve = TMP_ALLOC_LIMBS (primesieve_size (sieve_limit));
167 prime_limit = gmp_primesieve(sieve, sieve_limit);
168
169 /* TODO: Storing (prime - last_prime)/2 would allow this to go
170 up to the gap 304599508537+514=304599509051 .
171 With the current code our limit is 436273009+282=436273291 */
172 ASSERT (sieve_limit < 436273291);
173 /* THINK: Memory used by both sieve and primegap_tmp is kept
174 allocated, but they may overlap if primegap is filled from
175 larger down to smaller primes...
176 */
177
178 /* Needed to avoid assignment of read-only location */
179 primegap_tmp = TMP_ALLOC_TYPE (prime_limit, unsigned char);
180 primegap = primegap_tmp;
181
182 i = 0;
183 last_prime = 3;
184 /* THINK: should we get rid of sieve_limit and use (i < prime_limit)? */
185 for (mp_limb_t j = 4, *sp = sieve; j < sieve_limit; j += GMP_LIMB_BITS * 3)
186 for (mp_limb_t b = j, x = ~ *(sp++); x != 0; b += 3, x >>= 1)
187 if (x & 1)
188 {
189 mp_limb_t prime = b | 1;
190 primegap_tmp[i++] = prime - last_prime;
191 last_prime = prime;
192 }
193
194 /* Both primesieve and prime_limit ignore the first two primes. */
195 ASSERT(i == prime_limit);
196 }
197
198 if (nbits <= 32)
199 odds_in_composite_sieve = 336 / 2;
200 else if (nbits <= 64)
201 odds_in_composite_sieve = 1550 / 2;
202 else
203 /* Corresponds to a merit 14 prime_gap, which is rare. */
204 odds_in_composite_sieve = 5 * nbits;
205
206 /* composite[2*i] stores if p+2*i is a known composite */
207 composite = TMP_ALLOC_TYPE (odds_in_composite_sieve, char);
208
209 for (;;)
210 {
211 unsigned long difference;
212 unsigned long incr, prime;
213 int primetest;
214
215 memset (composite, 0, odds_in_composite_sieve);
216 prime = 3;
217 for (i = 0; i < prime_limit; i++)
218 {
219 /* Distance to next multiple of prime */
220 m = negative_mod_ui(p, prime);
221 /* Only care about odd multiplies of prime. */
222 if (m & 1)
223 m += prime;
224 m >>= 1;
225
226 /* Mark off any composites in sieve */
227 for (; m < odds_in_composite_sieve; m += prime)
228 composite[m] = 1;
229 prime += primegap[i];
230 }
231
232 difference = 0;
233 for (incr = 0; incr < odds_in_composite_sieve; difference += 2, incr += 1)
234 {
235 if (composite[incr])
236 continue;
237
238 increment_ui(p, p, difference);
239 difference = 0;
240
241 /* Miller-Rabin test */
242 primetest = mpz_millerrabin (p, 25);
243 if (primetest)
244 {
245 TMP_FREE;
246 return primetest;
247 }
248 }
249
250 /* Sieve next segment, very rare */
251 increment_ui(p, p, difference);
252 }
253 }
254
255 void
256 mpz_nextprime (mpz_ptr p, mpz_srcptr n)
257 {
258 /* Handle negative and small numbers */
259 if (mpz_cmp_ui (n, NP_SMALL_LIMIT) < 0)
260 {
261 ASSERT (NP_SMALL_LIMIT < UINT_MAX);
262 mpz_set_ui (p, findnext_small (SIZ (n) > 0 ? mpz_get_ui (n) : 1, +2));
263 return;
264 }
265
266 /* First odd greater than n */
267 mpz_add_ui (p, n, 1);
268
269 findnext(p, mpz_cdiv_ui, mpz_add_ui);
270 }
271
272 int
273 mpz_prevprime (mpz_ptr p, mpz_srcptr n)
274 {
275 /* Handle negative and small numbers */
276 if (mpz_cmp_ui (n, 2) <= 0)
277 return 0;
278
279 if (mpz_cmp_ui (n, NP_SMALL_LIMIT) < 0)
280 {
281 ASSERT (NP_SMALL_LIMIT < UINT_MAX);
282 mpz_set_ui (p, findnext_small (mpz_get_ui (n), -2));
283 return 2;
284 }
285
286 /* First odd less than n */
287 mpz_sub_ui (p, n, 2);
288
289 return findnext(p, mpz_tdiv_ui, mpz_sub_ui);
290 }
291