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