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 #ifdef __clang__ 15 #pragma clang diagnostic ignored "-Wtautological-constant-out-of-range-compare" 16 #endif 17 18 _LIBCPP_BEGIN_NAMESPACE_STD 19 20 namespace { 21 22 // handle all next_prime(i) for i in [1, 210), special case 0 23 const unsigned small_primes[] = 24 { 25 0, 26 2, 27 3, 28 5, 29 7, 30 11, 31 13, 32 17, 33 19, 34 23, 35 29, 36 31, 37 37, 38 41, 39 43, 40 47, 41 53, 42 59, 43 61, 44 67, 45 71, 46 73, 47 79, 48 83, 49 89, 50 97, 51 101, 52 103, 53 107, 54 109, 55 113, 56 127, 57 131, 58 137, 59 139, 60 149, 61 151, 62 157, 63 163, 64 167, 65 173, 66 179, 67 181, 68 191, 69 193, 70 197, 71 199, 72 211 73 }; 74 75 // potential primes = 210*k + indices[i], k >= 1 76 // these numbers are not divisible by 2, 3, 5 or 7 77 // (or any integer 2 <= j <= 10 for that matter). 78 const unsigned indices[] = 79 { 80 1, 81 11, 82 13, 83 17, 84 19, 85 23, 86 29, 87 31, 88 37, 89 41, 90 43, 91 47, 92 53, 93 59, 94 61, 95 67, 96 71, 97 73, 98 79, 99 83, 100 89, 101 97, 102 101, 103 103, 104 107, 105 109, 106 113, 107 121, 108 127, 109 131, 110 137, 111 139, 112 143, 113 149, 114 151, 115 157, 116 163, 117 167, 118 169, 119 173, 120 179, 121 181, 122 187, 123 191, 124 193, 125 197, 126 199, 127 209 128 }; 129 130 } 131 132 // Returns: If n == 0, returns 0. Else returns the lowest prime number that 133 // is greater than or equal to n. 134 // 135 // The algorithm creates a list of small primes, plus an open-ended list of 136 // potential primes. All prime numbers are potential prime numbers. However 137 // some potential prime numbers are not prime. In an ideal world, all potential 138 // prime numbers would be prime. Candidate prime numbers are chosen as the next 139 // highest potential prime. Then this number is tested for prime by dividing it 140 // by all potential prime numbers less than the sqrt of the candidate. 141 // 142 // This implementation defines potential primes as those numbers not divisible 143 // by 2, 3, 5, and 7. Other (common) implementations define potential primes 144 // as those not divisible by 2. A few other implementations define potential 145 // primes as those not divisible by 2 or 3. By raising the number of small 146 // primes which the potential prime is not divisible by, the set of potential 147 // primes more closely approximates the set of prime numbers. And thus there 148 // are fewer potential primes to search, and fewer potential primes to divide 149 // against. 150 151 template <size_t _Sz = sizeof(size_t)> 152 inline _LIBCPP_INLINE_VISIBILITY 153 typename enable_if<_Sz == 4, void>::type 154 __check_for_overflow(size_t N) 155 { 156 if (N > 0xFFFFFFFB) 157 __throw_overflow_error("__next_prime overflow"); 158 } 159 160 template <size_t _Sz = sizeof(size_t)> 161 inline _LIBCPP_INLINE_VISIBILITY 162 typename enable_if<_Sz == 8, void>::type 163 __check_for_overflow(size_t N) 164 { 165 if (N > 0xFFFFFFFFFFFFFFC5ull) 166 __throw_overflow_error("__next_prime overflow"); 167 } 168 169 size_t 170 __next_prime(size_t n) 171 { 172 const size_t L = 210; 173 const size_t N = sizeof(small_primes) / sizeof(small_primes[0]); 174 // If n is small enough, search in small_primes 175 if (n <= small_primes[N-1]) 176 return *std::lower_bound(small_primes, small_primes + N, n); 177 // Else n > largest small_primes 178 // Check for overflow 179 __check_for_overflow(n); 180 // Start searching list of potential primes: L * k0 + indices[in] 181 const size_t M = sizeof(indices) / sizeof(indices[0]); 182 // Select first potential prime >= n 183 // Known a-priori n >= L 184 size_t k0 = n / L; 185 size_t in = static_cast<size_t>(std::lower_bound(indices, indices + M, n - k0 * L) 186 - indices); 187 n = L * k0 + indices[in]; 188 while (true) 189 { 190 // Divide n by all primes or potential primes (i) until: 191 // 1. The division is even, so try next potential prime. 192 // 2. The i > sqrt(n), in which case n is prime. 193 // It is known a-priori that n is not divisible by 2, 3, 5 or 7, 194 // so don't test those (j == 5 -> divide by 11 first). And the 195 // potential primes start with 211, so don't test against the last 196 // small prime. 197 for (size_t j = 5; j < N - 1; ++j) 198 { 199 const std::size_t p = small_primes[j]; 200 const std::size_t q = n / p; 201 if (q < p) 202 return n; 203 if (n == q * p) 204 goto next; 205 } 206 // n wasn't divisible by small primes, try potential primes 207 { 208 size_t i = 211; 209 while (true) 210 { 211 std::size_t q = n / i; 212 if (q < i) 213 return n; 214 if (n == q * i) 215 break; 216 217 i += 10; 218 q = n / i; 219 if (q < i) 220 return n; 221 if (n == q * i) 222 break; 223 224 i += 2; 225 q = n / i; 226 if (q < i) 227 return n; 228 if (n == q * i) 229 break; 230 231 i += 4; 232 q = n / i; 233 if (q < i) 234 return n; 235 if (n == q * i) 236 break; 237 238 i += 2; 239 q = n / i; 240 if (q < i) 241 return n; 242 if (n == q * i) 243 break; 244 245 i += 4; 246 q = n / i; 247 if (q < i) 248 return n; 249 if (n == q * i) 250 break; 251 252 i += 6; 253 q = n / i; 254 if (q < i) 255 return n; 256 if (n == q * i) 257 break; 258 259 i += 2; 260 q = n / i; 261 if (q < i) 262 return n; 263 if (n == q * i) 264 break; 265 266 i += 6; 267 q = n / i; 268 if (q < i) 269 return n; 270 if (n == q * i) 271 break; 272 273 i += 4; 274 q = n / i; 275 if (q < i) 276 return n; 277 if (n == q * i) 278 break; 279 280 i += 2; 281 q = n / i; 282 if (q < i) 283 return n; 284 if (n == q * i) 285 break; 286 287 i += 4; 288 q = n / i; 289 if (q < i) 290 return n; 291 if (n == q * i) 292 break; 293 294 i += 6; 295 q = n / i; 296 if (q < i) 297 return n; 298 if (n == q * i) 299 break; 300 301 i += 6; 302 q = n / i; 303 if (q < i) 304 return n; 305 if (n == q * i) 306 break; 307 308 i += 2; 309 q = n / i; 310 if (q < i) 311 return n; 312 if (n == q * i) 313 break; 314 315 i += 6; 316 q = n / i; 317 if (q < i) 318 return n; 319 if (n == q * i) 320 break; 321 322 i += 4; 323 q = n / i; 324 if (q < i) 325 return n; 326 if (n == q * i) 327 break; 328 329 i += 2; 330 q = n / i; 331 if (q < i) 332 return n; 333 if (n == q * i) 334 break; 335 336 i += 6; 337 q = n / i; 338 if (q < i) 339 return n; 340 if (n == q * i) 341 break; 342 343 i += 4; 344 q = n / i; 345 if (q < i) 346 return n; 347 if (n == q * i) 348 break; 349 350 i += 6; 351 q = n / i; 352 if (q < i) 353 return n; 354 if (n == q * i) 355 break; 356 357 i += 8; 358 q = n / i; 359 if (q < i) 360 return n; 361 if (n == q * i) 362 break; 363 364 i += 4; 365 q = n / i; 366 if (q < i) 367 return n; 368 if (n == q * i) 369 break; 370 371 i += 2; 372 q = n / i; 373 if (q < i) 374 return n; 375 if (n == q * i) 376 break; 377 378 i += 4; 379 q = n / i; 380 if (q < i) 381 return n; 382 if (n == q * i) 383 break; 384 385 i += 2; 386 q = n / i; 387 if (q < i) 388 return n; 389 if (n == q * i) 390 break; 391 392 i += 4; 393 q = n / i; 394 if (q < i) 395 return n; 396 if (n == q * i) 397 break; 398 399 i += 8; 400 q = n / i; 401 if (q < i) 402 return n; 403 if (n == q * i) 404 break; 405 406 i += 6; 407 q = n / i; 408 if (q < i) 409 return n; 410 if (n == q * i) 411 break; 412 413 i += 4; 414 q = n / i; 415 if (q < i) 416 return n; 417 if (n == q * i) 418 break; 419 420 i += 6; 421 q = n / i; 422 if (q < i) 423 return n; 424 if (n == q * i) 425 break; 426 427 i += 2; 428 q = n / i; 429 if (q < i) 430 return n; 431 if (n == q * i) 432 break; 433 434 i += 4; 435 q = n / i; 436 if (q < i) 437 return n; 438 if (n == q * i) 439 break; 440 441 i += 6; 442 q = n / i; 443 if (q < i) 444 return n; 445 if (n == q * i) 446 break; 447 448 i += 2; 449 q = n / i; 450 if (q < i) 451 return n; 452 if (n == q * i) 453 break; 454 455 i += 6; 456 q = n / i; 457 if (q < i) 458 return n; 459 if (n == q * i) 460 break; 461 462 i += 6; 463 q = n / i; 464 if (q < i) 465 return n; 466 if (n == q * i) 467 break; 468 469 i += 4; 470 q = n / i; 471 if (q < i) 472 return n; 473 if (n == q * i) 474 break; 475 476 i += 2; 477 q = n / i; 478 if (q < i) 479 return n; 480 if (n == q * i) 481 break; 482 483 i += 4; 484 q = n / i; 485 if (q < i) 486 return n; 487 if (n == q * i) 488 break; 489 490 i += 6; 491 q = n / i; 492 if (q < i) 493 return n; 494 if (n == q * i) 495 break; 496 497 i += 2; 498 q = n / i; 499 if (q < i) 500 return n; 501 if (n == q * i) 502 break; 503 504 i += 6; 505 q = n / i; 506 if (q < i) 507 return n; 508 if (n == q * i) 509 break; 510 511 i += 4; 512 q = n / i; 513 if (q < i) 514 return n; 515 if (n == q * i) 516 break; 517 518 i += 2; 519 q = n / i; 520 if (q < i) 521 return n; 522 if (n == q * i) 523 break; 524 525 i += 4; 526 q = n / i; 527 if (q < i) 528 return n; 529 if (n == q * i) 530 break; 531 532 i += 2; 533 q = n / i; 534 if (q < i) 535 return n; 536 if (n == q * i) 537 break; 538 539 i += 10; 540 q = n / i; 541 if (q < i) 542 return n; 543 if (n == q * i) 544 break; 545 546 // This will loop i to the next "plane" of potential primes 547 i += 2; 548 } 549 } 550 next: 551 // n is not prime. Increment n to next potential prime. 552 if (++in == M) 553 { 554 ++k0; 555 in = 0; 556 } 557 n = L * k0 + indices[in]; 558 } 559 } 560 561 _LIBCPP_END_NAMESPACE_STD 562