1//===----------------------------------------------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#include <__hash_table>
10#include <algorithm>
11#include <stdexcept>
12
13_LIBCPP_CLANG_DIAGNOSTIC_IGNORED("-Wtautological-constant-out-of-range-compare")
14
15_LIBCPP_BEGIN_NAMESPACE_STD
16
17namespace {
18
19// handle all next_prime(i) for i in [1, 210), special case 0
20const unsigned small_primes[] = {
21 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
22 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
23 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211};
24
25// potential primes = 210*k + indices[i], k >= 1
26// these numbers are not divisible by 2, 3, 5 or 7
27// (or any integer 2 <= j <= 10 for that matter).
28const unsigned indices[] = {
29 1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
30 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139,
31 143, 149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209};
32
33} // namespace
34
35// Returns: If n == 0, returns 0. Else returns the lowest prime number that
36// is greater than or equal to n.
37//
38// The algorithm creates a list of small primes, plus an open-ended list of
39// potential primes. All prime numbers are potential prime numbers. However
40// some potential prime numbers are not prime. In an ideal world, all potential
41// prime numbers would be prime. Candidate prime numbers are chosen as the next
42// highest potential prime. Then this number is tested for prime by dividing it
43// by all potential prime numbers less than the sqrt of the candidate.
44//
45// This implementation defines potential primes as those numbers not divisible
46// by 2, 3, 5, and 7. Other (common) implementations define potential primes
47// as those not divisible by 2. A few other implementations define potential
48// primes as those not divisible by 2 or 3. By raising the number of small
49// primes which the potential prime is not divisible by, the set of potential
50// primes more closely approximates the set of prime numbers. And thus there
51// are fewer potential primes to search, and fewer potential primes to divide
52// against.
53
54inline void __check_for_overflow(size_t N) {
55 if constexpr (sizeof(size_t) == 4) {
56 if (N > 0xFFFFFFFB)
57 std::__throw_overflow_error(msg: "__next_prime overflow");
58 } else {
59 if (N > 0xFFFFFFFFFFFFFFC5ull)
60 std::__throw_overflow_error(msg: "__next_prime overflow");
61 }
62}
63
64size_t __next_prime(size_t n) {
65 const size_t L = 210;
66 const size_t N = sizeof(small_primes) / sizeof(small_primes[0]);
67 // If n is small enough, search in small_primes
68 if (n <= small_primes[N - 1])
69 return *std::lower_bound(first: small_primes, last: small_primes + N, value: n);
70 // Else n > largest small_primes
71 // Check for overflow
72 __check_for_overflow(N: n);
73 // Start searching list of potential primes: L * k0 + indices[in]
74 const size_t M = sizeof(indices) / sizeof(indices[0]);
75 // Select first potential prime >= n
76 // Known a-priori n >= L
77 size_t k0 = n / L;
78 size_t in = static_cast<size_t>(std::lower_bound(first: indices, last: indices + M, value: n - k0 * L) - indices);
79 n = L * k0 + indices[in];
80 while (true) {
81 // Divide n by all primes or potential primes (i) until:
82 // 1. The division is even, so try next potential prime.
83 // 2. The i > sqrt(n), in which case n is prime.
84 // It is known a-priori that n is not divisible by 2, 3, 5 or 7,
85 // so don't test those (j == 5 -> divide by 11 first). And the
86 // potential primes start with 211, so don't test against the last
87 // small prime.
88 for (size_t j = 5; j < N - 1; ++j) {
89 const std::size_t p = small_primes[j];
90 const std::size_t q = n / p;
91 if (q < p)
92 return n;
93 if (n == q * p)
94 goto next;
95 }
96 // n wasn't divisible by small primes, try potential primes
97 {
98 size_t i = 211;
99 while (true) {
100 std::size_t q = n / i;
101 if (q < i)
102 return n;
103 if (n == q * i)
104 break;
105
106 i += 10;
107 q = n / i;
108 if (q < i)
109 return n;
110 if (n == q * i)
111 break;
112
113 i += 2;
114 q = n / i;
115 if (q < i)
116 return n;
117 if (n == q * i)
118 break;
119
120 i += 4;
121 q = n / i;
122 if (q < i)
123 return n;
124 if (n == q * i)
125 break;
126
127 i += 2;
128 q = n / i;
129 if (q < i)
130 return n;
131 if (n == q * i)
132 break;
133
134 i += 4;
135 q = n / i;
136 if (q < i)
137 return n;
138 if (n == q * i)
139 break;
140
141 i += 6;
142 q = n / i;
143 if (q < i)
144 return n;
145 if (n == q * i)
146 break;
147
148 i += 2;
149 q = n / i;
150 if (q < i)
151 return n;
152 if (n == q * i)
153 break;
154
155 i += 6;
156 q = n / i;
157 if (q < i)
158 return n;
159 if (n == q * i)
160 break;
161
162 i += 4;
163 q = n / i;
164 if (q < i)
165 return n;
166 if (n == q * i)
167 break;
168
169 i += 2;
170 q = n / i;
171 if (q < i)
172 return n;
173 if (n == q * i)
174 break;
175
176 i += 4;
177 q = n / i;
178 if (q < i)
179 return n;
180 if (n == q * i)
181 break;
182
183 i += 6;
184 q = n / i;
185 if (q < i)
186 return n;
187 if (n == q * i)
188 break;
189
190 i += 6;
191 q = n / i;
192 if (q < i)
193 return n;
194 if (n == q * i)
195 break;
196
197 i += 2;
198 q = n / i;
199 if (q < i)
200 return n;
201 if (n == q * i)
202 break;
203
204 i += 6;
205 q = n / i;
206 if (q < i)
207 return n;
208 if (n == q * i)
209 break;
210
211 i += 4;
212 q = n / i;
213 if (q < i)
214 return n;
215 if (n == q * i)
216 break;
217
218 i += 2;
219 q = n / i;
220 if (q < i)
221 return n;
222 if (n == q * i)
223 break;
224
225 i += 6;
226 q = n / i;
227 if (q < i)
228 return n;
229 if (n == q * i)
230 break;
231
232 i += 4;
233 q = n / i;
234 if (q < i)
235 return n;
236 if (n == q * i)
237 break;
238
239 i += 6;
240 q = n / i;
241 if (q < i)
242 return n;
243 if (n == q * i)
244 break;
245
246 i += 8;
247 q = n / i;
248 if (q < i)
249 return n;
250 if (n == q * i)
251 break;
252
253 i += 4;
254 q = n / i;
255 if (q < i)
256 return n;
257 if (n == q * i)
258 break;
259
260 i += 2;
261 q = n / i;
262 if (q < i)
263 return n;
264 if (n == q * i)
265 break;
266
267 i += 4;
268 q = n / i;
269 if (q < i)
270 return n;
271 if (n == q * i)
272 break;
273
274 i += 2;
275 q = n / i;
276 if (q < i)
277 return n;
278 if (n == q * i)
279 break;
280
281 i += 4;
282 q = n / i;
283 if (q < i)
284 return n;
285 if (n == q * i)
286 break;
287
288 i += 8;
289 q = n / i;
290 if (q < i)
291 return n;
292 if (n == q * i)
293 break;
294
295 i += 6;
296 q = n / i;
297 if (q < i)
298 return n;
299 if (n == q * i)
300 break;
301
302 i += 4;
303 q = n / i;
304 if (q < i)
305 return n;
306 if (n == q * i)
307 break;
308
309 i += 6;
310 q = n / i;
311 if (q < i)
312 return n;
313 if (n == q * i)
314 break;
315
316 i += 2;
317 q = n / i;
318 if (q < i)
319 return n;
320 if (n == q * i)
321 break;
322
323 i += 4;
324 q = n / i;
325 if (q < i)
326 return n;
327 if (n == q * i)
328 break;
329
330 i += 6;
331 q = n / i;
332 if (q < i)
333 return n;
334 if (n == q * i)
335 break;
336
337 i += 2;
338 q = n / i;
339 if (q < i)
340 return n;
341 if (n == q * i)
342 break;
343
344 i += 6;
345 q = n / i;
346 if (q < i)
347 return n;
348 if (n == q * i)
349 break;
350
351 i += 6;
352 q = n / i;
353 if (q < i)
354 return n;
355 if (n == q * i)
356 break;
357
358 i += 4;
359 q = n / i;
360 if (q < i)
361 return n;
362 if (n == q * i)
363 break;
364
365 i += 2;
366 q = n / i;
367 if (q < i)
368 return n;
369 if (n == q * i)
370 break;
371
372 i += 4;
373 q = n / i;
374 if (q < i)
375 return n;
376 if (n == q * i)
377 break;
378
379 i += 6;
380 q = n / i;
381 if (q < i)
382 return n;
383 if (n == q * i)
384 break;
385
386 i += 2;
387 q = n / i;
388 if (q < i)
389 return n;
390 if (n == q * i)
391 break;
392
393 i += 6;
394 q = n / i;
395 if (q < i)
396 return n;
397 if (n == q * i)
398 break;
399
400 i += 4;
401 q = n / i;
402 if (q < i)
403 return n;
404 if (n == q * i)
405 break;
406
407 i += 2;
408 q = n / i;
409 if (q < i)
410 return n;
411 if (n == q * i)
412 break;
413
414 i += 4;
415 q = n / i;
416 if (q < i)
417 return n;
418 if (n == q * i)
419 break;
420
421 i += 2;
422 q = n / i;
423 if (q < i)
424 return n;
425 if (n == q * i)
426 break;
427
428 i += 10;
429 q = n / i;
430 if (q < i)
431 return n;
432 if (n == q * i)
433 break;
434
435 // This will loop i to the next "plane" of potential primes
436 i += 2;
437 }
438 }
439 next:
440 // n is not prime. Increment n to next potential prime.
441 if (++in == M) {
442 ++k0;
443 in = 0;
444 }
445 n = L * k0 + indices[in];
446 }
447}
448
449_LIBCPP_END_NAMESPACE_STD
450