mirror of
https://github.com/jart/cosmopolitan.git
synced 2025-01-31 11:37:35 +00:00
452 lines
10 KiB
C++
452 lines
10 KiB
C++
//===----------------------------------------------------------------------===//
|
|
//
|
|
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
|
|
// See https://llvm.org/LICENSE.txt for license information.
|
|
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
|
|
//
|
|
//===----------------------------------------------------------------------===//
|
|
|
|
#include <__hash_table>
|
|
#include <algorithm>
|
|
#include <stdexcept>
|
|
#include <type_traits>
|
|
|
|
_LIBCPP_CLANG_DIAGNOSTIC_IGNORED("-Wtautological-constant-out-of-range-compare")
|
|
|
|
_LIBCPP_BEGIN_NAMESPACE_STD
|
|
|
|
namespace {
|
|
|
|
// handle all next_prime(i) for i in [1, 210), special case 0
|
|
const unsigned small_primes[] = {
|
|
0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
|
|
53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
|
|
131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211};
|
|
|
|
// potential primes = 210*k + indices[i], k >= 1
|
|
// these numbers are not divisible by 2, 3, 5 or 7
|
|
// (or any integer 2 <= j <= 10 for that matter).
|
|
const unsigned indices[] = {
|
|
1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
|
|
71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139,
|
|
143, 149, 151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209};
|
|
|
|
} // namespace
|
|
|
|
// Returns: If n == 0, returns 0. Else returns the lowest prime number that
|
|
// is greater than or equal to n.
|
|
//
|
|
// The algorithm creates a list of small primes, plus an open-ended list of
|
|
// potential primes. All prime numbers are potential prime numbers. However
|
|
// some potential prime numbers are not prime. In an ideal world, all potential
|
|
// prime numbers would be prime. Candidate prime numbers are chosen as the next
|
|
// highest potential prime. Then this number is tested for prime by dividing it
|
|
// by all potential prime numbers less than the sqrt of the candidate.
|
|
//
|
|
// This implementation defines potential primes as those numbers not divisible
|
|
// by 2, 3, 5, and 7. Other (common) implementations define potential primes
|
|
// as those not divisible by 2. A few other implementations define potential
|
|
// primes as those not divisible by 2 or 3. By raising the number of small
|
|
// primes which the potential prime is not divisible by, the set of potential
|
|
// primes more closely approximates the set of prime numbers. And thus there
|
|
// are fewer potential primes to search, and fewer potential primes to divide
|
|
// against.
|
|
|
|
template <size_t _Sz = sizeof(size_t)>
|
|
inline _LIBCPP_HIDE_FROM_ABI typename enable_if<_Sz == 4, void>::type __check_for_overflow(size_t N) {
|
|
if (N > 0xFFFFFFFB)
|
|
__throw_overflow_error("__next_prime overflow");
|
|
}
|
|
|
|
template <size_t _Sz = sizeof(size_t)>
|
|
inline _LIBCPP_HIDE_FROM_ABI typename enable_if<_Sz == 8, void>::type __check_for_overflow(size_t N) {
|
|
if (N > 0xFFFFFFFFFFFFFFC5ull)
|
|
__throw_overflow_error("__next_prime overflow");
|
|
}
|
|
|
|
size_t __next_prime(size_t n) {
|
|
const size_t L = 210;
|
|
const size_t N = sizeof(small_primes) / sizeof(small_primes[0]);
|
|
// If n is small enough, search in small_primes
|
|
if (n <= small_primes[N - 1])
|
|
return *std::lower_bound(small_primes, small_primes + N, n);
|
|
// Else n > largest small_primes
|
|
// Check for overflow
|
|
__check_for_overflow(n);
|
|
// Start searching list of potential primes: L * k0 + indices[in]
|
|
const size_t M = sizeof(indices) / sizeof(indices[0]);
|
|
// Select first potential prime >= n
|
|
// Known a-priori n >= L
|
|
size_t k0 = n / L;
|
|
size_t in = static_cast<size_t>(std::lower_bound(indices, indices + M, n - k0 * L) - indices);
|
|
n = L * k0 + indices[in];
|
|
while (true) {
|
|
// Divide n by all primes or potential primes (i) until:
|
|
// 1. The division is even, so try next potential prime.
|
|
// 2. The i > sqrt(n), in which case n is prime.
|
|
// It is known a-priori that n is not divisible by 2, 3, 5 or 7,
|
|
// so don't test those (j == 5 -> divide by 11 first). And the
|
|
// potential primes start with 211, so don't test against the last
|
|
// small prime.
|
|
for (size_t j = 5; j < N - 1; ++j) {
|
|
const std::size_t p = small_primes[j];
|
|
const std::size_t q = n / p;
|
|
if (q < p)
|
|
return n;
|
|
if (n == q * p)
|
|
goto next;
|
|
}
|
|
// n wasn't divisible by small primes, try potential primes
|
|
{
|
|
size_t i = 211;
|
|
while (true) {
|
|
std::size_t q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 10;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 8;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 8;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 6;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 4;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 2;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
i += 10;
|
|
q = n / i;
|
|
if (q < i)
|
|
return n;
|
|
if (n == q * i)
|
|
break;
|
|
|
|
// This will loop i to the next "plane" of potential primes
|
|
i += 2;
|
|
}
|
|
}
|
|
next:
|
|
// n is not prime. Increment n to next potential prime.
|
|
if (++in == M) {
|
|
++k0;
|
|
in = 0;
|
|
}
|
|
n = L * k0 + indices[in];
|
|
}
|
|
}
|
|
|
|
_LIBCPP_END_NAMESPACE_STD
|