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