1193323Sed//===----------------------------------------------------------------------===// 2193323Sed// 3193323Sed// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4193323Sed// See https://llvm.org/LICENSE.txt for license information. 5193323Sed// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6193323Sed// 7193323Sed//===----------------------------------------------------------------------===// 8193323Sed 9193323Sed#include <__hash_table> 10193323Sed#include <algorithm> 11193323Sed#include <stdexcept> 12193323Sed#include <type_traits> 13193323Sed 14193323Sed_LIBCPP_CLANG_DIAGNOSTIC_IGNORED("-Wtautological-constant-out-of-range-compare") 15193323Sed 16193323Sed_LIBCPP_BEGIN_NAMESPACE_STD 17193323Sed 18193323Sednamespace { 19193323Sed 20193323Sed// handle all next_prime(i) for i in [1, 210), special case 0 21199989Srdivackyconst unsigned small_primes[] = 22193323Sed{ 23198090Srdivacky 0, 24193323Sed 2, 25193323Sed 3, 26193323Sed 5, 27193323Sed 7, 28193323Sed 11, 29193323Sed 13, 30193323Sed 17, 31193323Sed 19, 32193323Sed 23, 33239462Sdim 29, 34239462Sdim 31, 35239462Sdim 37, 36239462Sdim 41, 37239462Sdim 43, 38239462Sdim 47, 39239462Sdim 53, 40239462Sdim 59, 41239462Sdim 61, 42239462Sdim 67, 43198090Srdivacky 71, 44198090Srdivacky 73, 45198090Srdivacky 79, 46198090Srdivacky 83, 47198090Srdivacky 89, 48198090Srdivacky 97, 49198090Srdivacky 101, 50193323Sed 103, 51193323Sed 107, 52193323Sed 109, 53193323Sed 113, 54193323Sed 127, 55193323Sed 131, 56193323Sed 137, 57193323Sed 139, 58193323Sed 149, 59193323Sed 151, 60193323Sed 157, 61193323Sed 163, 62239462Sdim 167, 63193323Sed 173, 64193323Sed 179, 65193323Sed 181, 66193323Sed 191, 67193323Sed 193, 68193323Sed 197, 69193323Sed 199, 70193323Sed 211 71193323Sed}; 72193323Sed 73193323Sed// potential primes = 210*k + indices[i], k >= 1 74193323Sed// these numbers are not divisible by 2, 3, 5 or 7 75193323Sed// (or any integer 2 <= j <= 10 for that matter). 76193323Sedconst unsigned indices[] = 77193323Sed{ 78193323Sed 1, 79193323Sed 11, 80193323Sed 13, 81193323Sed 17, 82193323Sed 19, 83193323Sed 23, 84193323Sed 29, 85193323Sed 31, 86193323Sed 37, 87193323Sed 41, 88193323Sed 43, 89193323Sed 47, 90193323Sed 53, 91193323Sed 59, 92193323Sed 61, 93193323Sed 67, 94193323Sed 71, 95193323Sed 73, 96193323Sed 79, 97193323Sed 83, 98193323Sed 89, 99193323Sed 97, 100193323Sed 101, 101193323Sed 103, 102193323Sed 107, 103193323Sed 109, 104193323Sed 113, 105193323Sed 121, 106193323Sed 127, 107193323Sed 131, 108193323Sed 137, 109193323Sed 139, 110193323Sed 143, 111193323Sed 149, 112193323Sed 151, 113193323Sed 157, 114193323Sed 163, 115193323Sed 167, 116193323Sed 169, 117193323Sed 173, 118193323Sed 179, 119193323Sed 181, 120193323Sed 187, 121193323Sed 191, 122193323Sed 193, 123193323Sed 197, 124193323Sed 199, 125193323Sed 209 126193323Sed}; 127193323Sed 128193323Sed} 129193323Sed 130193323Sed// Returns: If n == 0, returns 0. Else returns the lowest prime number that 131193323Sed// is greater than or equal to n. 132193323Sed// 133193323Sed// The algorithm creates a list of small primes, plus an open-ended list of 134193323Sed// potential primes. All prime numbers are potential prime numbers. However 135193323Sed// some potential prime numbers are not prime. In an ideal world, all potential 136193323Sed// prime numbers would be prime. Candidate prime numbers are chosen as the next 137193323Sed// highest potential prime. Then this number is tested for prime by dividing it 138193323Sed// by all potential prime numbers less than the sqrt of the candidate. 139193323Sed// 140193323Sed// This implementation defines potential primes as those numbers not divisible 141193323Sed// by 2, 3, 5, and 7. Other (common) implementations define potential primes 142193323Sed// as those not divisible by 2. A few other implementations define potential 143193323Sed// primes as those not divisible by 2 or 3. By raising the number of small 144193323Sed// primes which the potential prime is not divisible by, the set of potential 145193323Sed// primes more closely approximates the set of prime numbers. And thus there 146193323Sed// are fewer potential primes to search, and fewer potential primes to divide 147193323Sed// against. 148193323Sed 149193323Sedtemplate <size_t _Sz = sizeof(size_t)> 150193323Sedinline _LIBCPP_INLINE_VISIBILITY 151193323Sedtypename enable_if<_Sz == 4, void>::type 152193323Sed__check_for_overflow(size_t N) 153193323Sed{ 154193323Sed if (N > 0xFFFFFFFB) 155193323Sed __throw_overflow_error("__next_prime overflow"); 156193323Sed} 157193323Sed 158193323Sedtemplate <size_t _Sz = sizeof(size_t)> 159193323Sedinline _LIBCPP_INLINE_VISIBILITY 160193323Sedtypename enable_if<_Sz == 8, void>::type 161193323Sed__check_for_overflow(size_t N) 162193323Sed{ 163193323Sed if (N > 0xFFFFFFFFFFFFFFC5ull) 164193323Sed __throw_overflow_error("__next_prime overflow"); 165193323Sed} 166193323Sed 167193323Sedsize_t 168193323Sed__next_prime(size_t n) 169193323Sed{ 170193323Sed const size_t L = 210; 171193323Sed const size_t N = sizeof(small_primes) / sizeof(small_primes[0]); 172193323Sed // If n is small enough, search in small_primes 173193323Sed if (n <= small_primes[N-1]) 174193323Sed return *std::lower_bound(small_primes, small_primes + N, n); 175193323Sed // Else n > largest small_primes 176193323Sed // Check for overflow 177193323Sed __check_for_overflow(n); 178193323Sed // Start searching list of potential primes: L * k0 + indices[in] 179193323Sed const size_t M = sizeof(indices) / sizeof(indices[0]); 180193323Sed // Select first potential prime >= n 181193323Sed // Known a-priori n >= L 182193323Sed size_t k0 = n / L; 183193323Sed size_t in = static_cast<size_t>(std::lower_bound(indices, indices + M, n - k0 * L) 184193323Sed - indices); 185193323Sed n = L * k0 + indices[in]; 186193323Sed while (true) 187193323Sed { 188193323Sed // Divide n by all primes or potential primes (i) until: 189193323Sed // 1. The division is even, so try next potential prime. 190193323Sed // 2. The i > sqrt(n), in which case n is prime. 191193323Sed // It is known a-priori that n is not divisible by 2, 3, 5 or 7, 192193323Sed // so don't test those (j == 5 -> divide by 11 first). And the 193193323Sed // potential primes start with 211, so don't test against the last 194193323Sed // small prime. 195193323Sed for (size_t j = 5; j < N - 1; ++j) 196193323Sed { 197193323Sed const std::size_t p = small_primes[j]; 198193323Sed const std::size_t q = n / p; 199226633Sdim if (q < p) 200226633Sdim return n; 201226633Sdim if (n == q * p) 202226633Sdim goto next; 203193323Sed } 204226633Sdim // n wasn't divisible by small primes, try potential primes 205226633Sdim { 206193323Sed size_t i = 211; 207226633Sdim while (true) 208226633Sdim { 209226633Sdim std::size_t q = n / i; 210226633Sdim if (q < i) 211226633Sdim return n; 212226633Sdim if (n == q * i) 213226633Sdim break; 214193323Sed 215193323Sed i += 10; 216193323Sed q = n / i; 217193323Sed if (q < i) 218193323Sed return n; 219193323Sed if (n == q * i) 220263508Sdim break; 221263508Sdim 222263508Sdim i += 2; 223263508Sdim q = n / i; 224263508Sdim if (q < i) 225263508Sdim return n; 226263508Sdim if (n == q * i) 227263508Sdim break; 228263508Sdim 229263508Sdim i += 4; 230263508Sdim q = n / i; 231263508Sdim if (q < i) 232263508Sdim return n; 233263508Sdim if (n == q * i) 234263508Sdim break; 235263508Sdim 236193323Sed i += 2; 237193323Sed q = n / i; 238193323Sed if (q < i) 239193323Sed return n; 240193323Sed if (n == q * i) 241193323Sed break; 242193323Sed 243193323Sed i += 4; 244193323Sed q = n / i; 245193323Sed if (q < i) 246193323Sed return n; 247193323Sed if (n == q * i) 248193323Sed break; 249193323Sed 250212904Sdim i += 6; 251193323Sed q = n / i; 252193323Sed if (q < i) 253193323Sed return n; 254193323Sed if (n == q * i) 255193323Sed break; 256193323Sed 257239462Sdim i += 2; 258193323Sed q = n / i; 259193323Sed if (q < i) 260193323Sed return n; 261193323Sed if (n == q * i) 262193323Sed break; 263193323Sed 264193323Sed i += 6; 265249423Sdim q = n / i; 266193323Sed if (q < i) 267193323Sed return n; 268249423Sdim if (n == q * i) 269193323Sed break; 270193323Sed 271193323Sed i += 4; 272193323Sed q = n / i; 273193323Sed if (q < i) 274193323Sed return n; 275193323Sed if (n == q * i) 276193323Sed break; 277193323Sed 278193323Sed i += 2; 279193323Sed q = n / i; 280193323Sed if (q < i) 281193323Sed return n; 282193323Sed if (n == q * i) 283193323Sed break; 284193323Sed 285193323Sed i += 4; 286193323Sed q = n / i; 287193323Sed if (q < i) 288193323Sed return n; 289239462Sdim if (n == q * i) 290193323Sed break; 291193323Sed 292193323Sed i += 6; 293249423Sdim q = n / i; 294193323Sed if (q < i) 295193323Sed return n; 296263508Sdim if (n == q * i) 297263508Sdim break; 298263508Sdim 299263508Sdim i += 6; 300263508Sdim q = n / i; 301263508Sdim if (q < i) 302199481Srdivacky return n; 303199481Srdivacky if (n == q * i) 304263508Sdim break; 305263508Sdim 306199481Srdivacky i += 2; 307239462Sdim q = n / i; 308205218Srdivacky if (q < i) 309205218Srdivacky return n; 310205218Srdivacky if (n == q * i) 311205218Srdivacky break; 312205218Srdivacky 313205218Srdivacky i += 6; 314205218Srdivacky q = n / i; 315205218Srdivacky if (q < i) 316205218Srdivacky return n; 317205218Srdivacky if (n == q * i) 318205218Srdivacky break; 319205218Srdivacky 320205218Srdivacky i += 4; 321205218Srdivacky q = n / i; 322205218Srdivacky if (q < i) 323205218Srdivacky return n; 324205218Srdivacky if (n == q * i) 325205218Srdivacky break; 326205218Srdivacky 327205218Srdivacky i += 2; 328205218Srdivacky q = n / i; 329205218Srdivacky if (q < i) 330193323Sed return n; 331193323Sed if (n == q * i) 332193323Sed 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