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 17 namespace { 18 19 // handle all next_prime(i) for i in [1, 210), special case 0 20 const 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 >= 1 26 // these numbers are not divisible by 2, 3, 5 or 7 27 // (or any integer 2 <= j <= 10 for that matter). 28 const 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 } // namespace 34 35 // Returns: If n == 0, returns 0. Else returns the lowest prime number that 36 // is greater than or equal to n. 37 // 38 // The algorithm creates a list of small primes, plus an open-ended list of 39 // potential primes. All prime numbers are potential prime numbers. However 40 // some potential prime numbers are not prime. In an ideal world, all potential 41 // prime numbers would be prime. Candidate prime numbers are chosen as the next 42 // highest potential prime. Then this number is tested for prime by dividing it 43 // by all potential prime numbers less than the sqrt of the candidate. 44 // 45 // This implementation defines potential primes as those numbers not divisible 46 // by 2, 3, 5, and 7. Other (common) implementations define potential primes 47 // as those not divisible by 2. A few other implementations define potential 48 // primes as those not divisible by 2 or 3. By raising the number of small 49 // primes which the potential prime is not divisible by, the set of potential 50 // primes more closely approximates the set of prime numbers. And thus there 51 // are fewer potential primes to search, and fewer potential primes to divide 52 // against. 53 54 inline 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 64 size_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_primes 68 if (n <= small_primes[N - 1]) 69 return *std::lower_bound(small_primes, small_primes + N, n); 70 // Else n > largest small_primes 71 // Check for overflow 72 __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 >= n 76 // Known a-priori n >= L 77 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 the 86 // potential primes start with 211, so don't test against the last 87 // 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 primes 97 { 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 primes 436 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_STD 450