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 { 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). 76 const 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 149 template <size_t _Sz = sizeof(size_t)> 150 inline _LIBCPP_INLINE_VISIBILITY 151 typename 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 158 template <size_t _Sz = sizeof(size_t)> 159 inline _LIBCPP_INLINE_VISIBILITY 160 typename 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 167 size_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 } 548 next: 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