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