blob: 96cbd38f4dd84c708622ae4a146823d739edb262 [file] [log] [blame]
//===----------------------------------------------------------------------===//
//
// 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
};
}
// 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_INLINE_VISIBILITY
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_INLINE_VISIBILITY
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