450 lines · cpp
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-exception6//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_STD16 17namespace {18 19// handle all next_prime(i) for i in [1, 210), special case 020const 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 >= 126// these numbers are not divisible by 2, 3, 5 or 727// (or any integer 2 <= j <= 10 for that matter).28const 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} // namespace34 35// Returns: If n == 0, returns 0. Else returns the lowest prime number that36// is greater than or equal to n.37//38// The algorithm creates a list of small primes, plus an open-ended list of39// potential primes. All prime numbers are potential prime numbers. However40// some potential prime numbers are not prime. In an ideal world, all potential41// prime numbers would be prime. Candidate prime numbers are chosen as the next42// highest potential prime. Then this number is tested for prime by dividing it43// by all potential prime numbers less than the sqrt of the candidate.44//45// This implementation defines potential primes as those numbers not divisible46// by 2, 3, 5, and 7. Other (common) implementations define potential primes47// as those not divisible by 2. A few other implementations define potential48// primes as those not divisible by 2 or 3. By raising the number of small49// primes which the potential prime is not divisible by, the set of potential50// primes more closely approximates the set of prime numbers. And thus there51// are fewer potential primes to search, and fewer potential primes to divide52// against.53 54inline 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 64size_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_primes68 if (n <= small_primes[N - 1])69 return *std::lower_bound(small_primes, small_primes + N, n);70 // Else n > largest small_primes71 // Check for overflow72 __check_for_overflow(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 >= n76 // Known a-priori n >= L77 size_t k0 = n / L;78 size_t in = static_cast<size_t>(std::lower_bound(indices, indices + M, 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 the86 // potential primes start with 211, so don't test against the last87 // 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 primes97 {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 primes436 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_STD450