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)>
__check_for_overflow(size_t N)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)>
__check_for_overflow(size_t N)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
__next_prime(size_t n)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(small_primes, small_primes + N, 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(indices, indices + M, 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