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_LIBCPP_BEGIN_EXPLICIT_ABI_ANNOTATIONS
17
18namespace {
19
20// handle all next_prime(i) for i in [1, 210), special case 0
21const 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).
29const 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// These are the amount we increment by when checking for potential
35// primes in the loop in __next_prime.
36const uint8_t increments[] = {
37 0, 10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2,
38 4, 2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10,
39};
40
41} // namespace
42
43// Returns: If n == 0, returns 0. Else returns the lowest prime number that
44// is greater than or equal to n.
45//
46// The algorithm creates a list of small primes, plus an open-ended list of
47// potential primes. All prime numbers are potential prime numbers. However
48// some potential prime numbers are not prime. In an ideal world, all potential
49// prime numbers would be prime. Candidate prime numbers are chosen as the next
50// highest potential prime. Then this number is tested for prime by dividing it
51// by all potential prime numbers less than the sqrt of the candidate.
52//
53// This implementation defines potential primes as those numbers not divisible
54// by 2, 3, 5, and 7. Other (common) implementations define potential primes
55// as those not divisible by 2. A few other implementations define potential
56// primes as those not divisible by 2 or 3. By raising the number of small
57// primes which the potential prime is not divisible by, the set of potential
58// primes more closely approximates the set of prime numbers. And thus there
59// are fewer potential primes to search, and fewer potential primes to divide
60// against.
61
62inline void __check_for_overflow(size_t N) {
63 if constexpr (sizeof(size_t) == 4) {
64 if (N > 0xFFFFFFFB)
65 std::__throw_overflow_error(msg: "__next_prime overflow");
66 } else {
67 if (N > 0xFFFFFFFFFFFFFFC5ull)
68 std::__throw_overflow_error(msg: "__next_prime overflow");
69 }
70}
71
72size_t __next_prime(size_t n) {
73 const size_t L = 210;
74 const size_t N = sizeof(small_primes) / sizeof(small_primes[0]);
75 // If n is small enough, search in small_primes
76 if (n <= small_primes[N - 1])
77 return *std::lower_bound(first: small_primes, last: small_primes + N, value: n);
78 // Else n > largest small_primes
79 // Check for overflow
80 __check_for_overflow(N: n);
81 // Start searching list of potential primes: L * k0 + indices[in]
82 const size_t M = sizeof(indices) / sizeof(indices[0]);
83 // Select first potential prime >= n
84 // Known a-priori n >= L
85 size_t k0 = n / L;
86 size_t in = static_cast<size_t>(std::lower_bound(first: indices, last: indices + M, value: n - k0 * L) - indices);
87 n = L * k0 + indices[in];
88 while (true) {
89 // Divide n by all primes or potential primes (i) until:
90 // 1. The division is even, so try next potential prime.
91 // 2. The i > sqrt(n), in which case n is prime.
92 // It is known a-priori that n is not divisible by 2, 3, 5 or 7,
93 // so don't test those (j == 5 -> divide by 11 first). And the
94 // potential primes start with 211, so don't test against the last
95 // small prime.
96 for (size_t j = 5; j < N - 1; ++j) {
97 const std::size_t p = small_primes[j];
98 const std::size_t q = n / p;
99 if (q < p)
100 return n;
101 if (n == q * p)
102 goto next;
103 }
104 // n wasn't divisible by small primes, try potential primes
105 {
106 size_t i = 211;
107 while (true) {
108 for (auto inc : increments) {
109 i += inc;
110 std::size_t q = n / i;
111 if (q < i)
112 return n;
113 if (n == q * i)
114 goto next;
115 }
116
117 // This will loop i to the next "plane" of potential primes
118 i += 2;
119 }
120 }
121 next:
122 // n is not prime. Increment n to next potential prime.
123 if (++in == M) {
124 ++k0;
125 in = 0;
126 }
127 n = L * k0 + indices[in];
128 }
129}
130
131_LIBCPP_END_EXPLICIT_ABI_ANNOTATIONS
132_LIBCPP_END_NAMESPACE_STD
133