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 | |
17 | namespace { |
18 | |
19 | // handle all next_prime(i) for i in [1, 210), special case 0 |
20 | const 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). |
28 | const 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 | |
54 | inline void __check_for_overflow(size_t N) { |
55 | if constexpr (sizeof(size_t) == 4) { |
56 | if (N > 0xFFFFFFFB) |
57 | std::__throw_overflow_error("__next_prime overflow" ); |
58 | } else { |
59 | if (N > 0xFFFFFFFFFFFFFFC5ull) |
60 | std::__throw_overflow_error("__next_prime overflow" ); |
61 | } |
62 | } |
63 | |
64 | size_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, val: 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, val: 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 | |