mirror of
https://github.com/mpusz/mp-units.git
synced 2025-06-25 01:01:33 +02:00
Merge branch 'master' of github.com:mpusz/units
This commit is contained in:
Mateusz Pusz
@ -100,7 +100,7 @@ namespace mp_units::detail {
|
||||
return add_mod(
|
||||
// Transform into "negative space" to make the first parameter as small as possible;
|
||||
// then, transform back.
|
||||
n - mul_mod(n % a, num_batches, n),
|
||||
(n - mul_mod(n % a, num_batches, n)) % n,
|
||||
|
||||
// Handle the leftover product (which is guaranteed to fit in the integer type).
|
||||
(a * (b % batch_size)) % n,
|
||||
@ -183,28 +183,222 @@ struct NumberDecomposition {
|
||||
return false;
|
||||
}
|
||||
|
||||
[[nodiscard]] consteval bool is_prime_by_trial_division(std::uintmax_t n)
|
||||
// The Jacobi symbol, notated as `(a/n)`, is defined for odd positive `n` and any integer `a`, taking values
|
||||
// in the set `{-1, 0, 1}`. Besides being a completely multiplicative function (so that, for example, both
|
||||
// (a*b/n) = (a/n) * (b/n), and (a/n*m) = (a/n) * (a/m)), it obeys the following symmetry rules, which enable
|
||||
// its calculation:
|
||||
//
|
||||
// 1. (a/1) = 1, and (1/n) = 1, for all a and n.
|
||||
//
|
||||
// 2. (a/n) = 0 whenever a and n have a nontrivial common factor.
|
||||
//
|
||||
// 3. (a/n) = (b/n) whenever (a % n) = (b % n).
|
||||
//
|
||||
// 4. (2a/n) = (a/n) if n % 8 = 1 or 7, and -(a/n) if n % 8 = 3 or 5.
|
||||
//
|
||||
// 5. (a/n) = (n/a) * x if a and n are both odd, positive, and coprime. Here, x is 1 if either (a % 4) = 1
|
||||
// or (n % 4) = 1, and -1 otherwise.
|
||||
//
|
||||
// 6. (-1/n) = 1 if n % 4 = 1, and -1 if n % 4 = 3.
|
||||
[[nodiscard]] consteval int jacobi_symbol(std::int64_t raw_a, std::uint64_t n)
|
||||
{
|
||||
for (std::uintmax_t f = 2; f * f <= n; f += 1 + (f % 2)) {
|
||||
if (n % f == 0) {
|
||||
return false;
|
||||
}
|
||||
// Rule 1: n=1 case.
|
||||
if (n == 1u) {
|
||||
return 1;
|
||||
}
|
||||
return true;
|
||||
|
||||
// Starting conditions: transform `a` to strictly non-negative values, setting `result` to the sign that we
|
||||
// pick up (if any) from following these rules (i.e., rules 3 and 6).
|
||||
int result = ((raw_a >= 0) || (n % 4u == 1u)) ? 1 : -1;
|
||||
auto a = static_cast<std::uint64_t>(raw_a < 0 ? -raw_a : raw_a) % n;
|
||||
|
||||
while (a != 0u) {
|
||||
// Rule 4.
|
||||
const int sign_for_even = (n % 8u == 1u || n % 8u == 7u) ? 1 : -1;
|
||||
while (a % 2u == 0u) {
|
||||
a /= 2u;
|
||||
result *= sign_for_even;
|
||||
}
|
||||
|
||||
// Rule 1: a=1 case.
|
||||
if (a == 1u) {
|
||||
return result;
|
||||
}
|
||||
|
||||
// Rule 2.
|
||||
if (std::gcd(a, n) != 1u) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Note that at this point, we know that `a` and `n` are coprime, and are both odd and positive.
|
||||
// Therefore, we meet the preconditions for rule 5 (the "flip-and-reduce" rule).
|
||||
result *= (n % 4u == 1u || a % 4u == 1u) ? 1 : -1;
|
||||
const std::uint64_t new_a = n % a;
|
||||
n = a;
|
||||
a = new_a;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
// Return the first factor of n, as long as it is either k or n.
|
||||
//
|
||||
// Precondition: no integer smaller than k evenly divides n.
|
||||
[[nodiscard]] constexpr std::optional<std::uintmax_t> first_factor_maybe(std::uintmax_t n, std::uintmax_t k)
|
||||
[[nodiscard]] consteval bool is_perfect_square(std::uint64_t n)
|
||||
{
|
||||
if (n % k == 0) {
|
||||
return k;
|
||||
if (n < 2u) {
|
||||
return true;
|
||||
}
|
||||
if (k * k > n) {
|
||||
return n;
|
||||
|
||||
std::uint64_t prev = n / 2u;
|
||||
while (true) {
|
||||
const std::uint64_t curr = (prev + n / prev) / 2u;
|
||||
if (curr * curr == n) {
|
||||
return true;
|
||||
}
|
||||
if (curr >= prev) {
|
||||
return false;
|
||||
}
|
||||
prev = curr;
|
||||
}
|
||||
return std::nullopt;
|
||||
}
|
||||
|
||||
// The "D" parameter in the Strong Lucas Probable Prime test.
|
||||
//
|
||||
// Default construction produces the first value to try, according to Selfridge's parameter selection.
|
||||
// Calling `successor()` on this repeatedly will produce the sequence of values to try.
|
||||
struct LucasDParameter {
|
||||
std::uint64_t mag = 5u;
|
||||
bool pos = true;
|
||||
|
||||
friend constexpr int as_int(LucasDParameter p)
|
||||
{
|
||||
int D = static_cast<int>(p.mag);
|
||||
return p.pos ? D : -D;
|
||||
}
|
||||
friend constexpr LucasDParameter successor(LucasDParameter p) { return {.mag = p.mag + 2u, .pos = !p.pos}; }
|
||||
};
|
||||
|
||||
// The first `D` in the infinite sequence {5, -7, 9, -11, ...} whose Jacobi symbol is -1. This is the D we
|
||||
// want to use for the Strong Lucas Probable Prime test.
|
||||
//
|
||||
// Precondition: `n` is not a perfect square.
|
||||
[[nodiscard]] consteval LucasDParameter find_first_D_with_jacobi_symbol_neg_one(std::uint64_t n)
|
||||
{
|
||||
LucasDParameter D{};
|
||||
while (jacobi_symbol(as_int(D), n) != -1) {
|
||||
D = successor(D);
|
||||
}
|
||||
return D;
|
||||
}
|
||||
|
||||
// An element of the Lucas sequence, with parameters tuned for the Strong Lucas test.
|
||||
//
|
||||
// The default values give the first element (i.e., k=1) of the sequence.
|
||||
struct LucasSequenceElement {
|
||||
std::uint64_t u = 1u;
|
||||
std::uint64_t v = 1u;
|
||||
};
|
||||
|
||||
// Produce the Lucas element whose index is twice the input element's index.
|
||||
[[nodiscard]] consteval LucasSequenceElement double_strong_lucas_index(const LucasSequenceElement& element,
|
||||
std::uint64_t n, LucasDParameter D)
|
||||
{
|
||||
const auto& [u, v] = element;
|
||||
|
||||
std::uint64_t v_squared = mul_mod(v, v, n);
|
||||
std::uint64_t D_u_squared = mul_mod(D.mag % n, mul_mod(u, u, n), n);
|
||||
std::uint64_t new_v = D.pos ? add_mod(v_squared, D_u_squared, n) : sub_mod(v_squared, D_u_squared, n);
|
||||
new_v = half_mod_odd(new_v, n);
|
||||
|
||||
return {.u = mul_mod(u, v, n), .v = new_v};
|
||||
}
|
||||
|
||||
[[nodiscard]] consteval LucasSequenceElement increment_strong_lucas_index(const LucasSequenceElement& element,
|
||||
std::uint64_t n, LucasDParameter D)
|
||||
{
|
||||
const auto& [u, v] = element;
|
||||
|
||||
const auto new_u = half_mod_odd(add_mod(u, v, n), n);
|
||||
|
||||
const auto D_u = mul_mod(D.mag % n, u, n);
|
||||
auto new_v = D.pos ? add_mod(v, D_u, n) : sub_mod(v, D_u, n);
|
||||
new_v = half_mod_odd(new_v, n);
|
||||
|
||||
return {.u = new_u, .v = new_v};
|
||||
}
|
||||
|
||||
[[nodiscard]] consteval LucasSequenceElement find_strong_lucas_element(std::uint64_t i, std::uint64_t n,
|
||||
LucasDParameter D)
|
||||
{
|
||||
LucasSequenceElement element{};
|
||||
|
||||
bool bits[64] = {};
|
||||
std::size_t n_bits = 0u;
|
||||
while (i > 1u) {
|
||||
bits[n_bits++] = (i & 1u);
|
||||
i >>= 1u;
|
||||
}
|
||||
|
||||
for (auto j = n_bits; j > 0u; --j) {
|
||||
element = double_strong_lucas_index(element, n, D);
|
||||
if (bits[j - 1u]) {
|
||||
element = increment_strong_lucas_index(element, n, D);
|
||||
}
|
||||
}
|
||||
|
||||
return element;
|
||||
}
|
||||
|
||||
// Perform a Strong Lucas Probable Prime test on `n`.
|
||||
//
|
||||
// Precondition: (n >= 2).
|
||||
// Precondition: (n is odd).
|
||||
[[nodiscard]] consteval bool strong_lucas_probable_prime(std::uint64_t n)
|
||||
{
|
||||
MP_UNITS_EXPECTS_DEBUG(n >= 2u);
|
||||
MP_UNITS_EXPECTS_DEBUG(n % 2u == 1u);
|
||||
|
||||
if (is_perfect_square(n)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
const auto D = find_first_D_with_jacobi_symbol_neg_one(n);
|
||||
|
||||
const auto [s, d] = decompose(n + 1u);
|
||||
|
||||
auto element = find_strong_lucas_element(d, n, D);
|
||||
if (element.u == 0u) {
|
||||
return true;
|
||||
}
|
||||
|
||||
for (auto i = 0u; i < s; ++i) {
|
||||
if (element.v == 0u) {
|
||||
return true;
|
||||
}
|
||||
element = double_strong_lucas_index(element, n, D);
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
// The Baillie-PSW test is technically a "probable prime" test. However, it is known to be correct for all
|
||||
// 64-bit integers, and no counterexample of any size has ever been found. Thus, for this library's purposes,
|
||||
// we can treat it as deterministic... probably.
|
||||
[[nodiscard]] consteval bool baillie_psw_probable_prime(std::uint64_t n)
|
||||
{
|
||||
if (n < 2u) {
|
||||
return false;
|
||||
}
|
||||
if (n < 4u) {
|
||||
return true;
|
||||
}
|
||||
if (n % 2u == 0u) {
|
||||
return false;
|
||||
}
|
||||
|
||||
if (!miller_rabin_probable_prime(2u, n)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
return strong_lucas_probable_prime(n);
|
||||
}
|
||||
|
||||
template<std::size_t N>
|
||||
@ -214,42 +408,13 @@ template<std::size_t N>
|
||||
primes[0] = 2;
|
||||
for (std::size_t i = 1; i < N; ++i) {
|
||||
primes[i] = primes[i - 1] + 1;
|
||||
while (!is_prime_by_trial_division(primes[i])) {
|
||||
while (!baillie_psw_probable_prime(primes[i])) {
|
||||
++primes[i];
|
||||
}
|
||||
}
|
||||
return primes;
|
||||
}
|
||||
|
||||
template<std::size_t N, typename Callable>
|
||||
consteval void call_for_coprimes_up_to(std::uintmax_t n, const std::array<std::uintmax_t, N>& basis, Callable call)
|
||||
{
|
||||
for (std::uintmax_t i = 0u; i < n; ++i) {
|
||||
if (std::apply([&i](auto... primes) { return ((i % primes != 0) && ...); }, basis)) {
|
||||
call(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<std::size_t N>
|
||||
[[nodiscard]] consteval std::size_t num_coprimes_up_to(std::uintmax_t n, const std::array<std::uintmax_t, N>& basis)
|
||||
{
|
||||
std::size_t count = 0u;
|
||||
call_for_coprimes_up_to(n, basis, [&count](auto) { ++count; });
|
||||
return count;
|
||||
}
|
||||
|
||||
template<std::size_t ResultSize, std::size_t N>
|
||||
[[nodiscard]] consteval auto coprimes_up_to(std::size_t n, const std::array<std::uintmax_t, N>& basis)
|
||||
{
|
||||
std::array<std::uintmax_t, ResultSize> coprimes{};
|
||||
std::size_t i = 0u;
|
||||
|
||||
call_for_coprimes_up_to(n, basis, [&coprimes, &i](std::uintmax_t cp) { coprimes[i++] = cp; });
|
||||
|
||||
return coprimes;
|
||||
}
|
||||
|
||||
template<std::size_t N>
|
||||
[[nodiscard]] consteval std::uintmax_t product(const std::array<std::uintmax_t, N>& values)
|
||||
{
|
||||
@ -276,48 +441,34 @@ constexpr auto get_first_of(const Rng& rng, UnaryFunction f)
|
||||
return get_first_of(begin(rng), end(rng), f);
|
||||
}
|
||||
|
||||
// A configurable instantiation of the "wheel factorization" algorithm [1] for prime numbers.
|
||||
//
|
||||
// Instantiate with N to use a "basis" of the first N prime numbers. Higher values of N use fewer trial divisions, at
|
||||
// the cost of additional space. The amount of space consumed is roughly the total number of numbers that are a) less
|
||||
// than the _product_ of the basis elements (first N primes), and b) coprime with every element of the basis. This
|
||||
// means it grows rapidly with N. Consider this approximate chart:
|
||||
//
|
||||
// N | Num coprimes | Trial divisions needed
|
||||
// --+--------------+-----------------------
|
||||
// 1 | 1 | 50.0 %
|
||||
// 2 | 2 | 33.3 %
|
||||
// 3 | 8 | 26.7 %
|
||||
// 4 | 48 | 22.9 %
|
||||
// 5 | 480 | 20.8 %
|
||||
//
|
||||
// Note the diminishing returns, and the rapidly escalating costs. Consider this behaviour when choosing the value of N
|
||||
// most appropriate for your needs.
|
||||
//
|
||||
// [1] https://en.wikipedia.org/wiki/Wheel_factorization
|
||||
template<std::size_t BasisSize>
|
||||
struct wheel_factorizer {
|
||||
static constexpr auto basis = first_n_primes<BasisSize>();
|
||||
static constexpr auto wheel_size = product(basis);
|
||||
static constexpr auto coprimes_in_first_wheel =
|
||||
coprimes_up_to<num_coprimes_up_to(wheel_size, basis)>(wheel_size, basis);
|
||||
[[nodiscard]] consteval std::uintmax_t find_first_factor(std::uintmax_t n)
|
||||
{
|
||||
constexpr auto first_100_primes = first_n_primes<100>();
|
||||
|
||||
[[nodiscard]] static consteval std::uintmax_t find_first_factor(std::uintmax_t n)
|
||||
{
|
||||
if (const auto k = detail::get_first_of(basis, [&](auto p) { return first_factor_maybe(n, p); })) return *k;
|
||||
for (const auto& p : first_100_primes) {
|
||||
if (n % p == 0u) {
|
||||
return p;
|
||||
}
|
||||
if (p * p > n) {
|
||||
return n;
|
||||
}
|
||||
}
|
||||
|
||||
if (const auto k = detail::get_first_of(std::next(begin(coprimes_in_first_wheel)), end(coprimes_in_first_wheel),
|
||||
[&](auto p) { return first_factor_maybe(n, p); }))
|
||||
return *k;
|
||||
|
||||
for (std::size_t wheel = wheel_size; wheel < n; wheel += wheel_size)
|
||||
if (const auto k =
|
||||
detail::get_first_of(coprimes_in_first_wheel, [&](auto p) { return first_factor_maybe(n, wheel + p); }))
|
||||
return *k;
|
||||
// If we got this far, and we haven't found a factor, and haven't terminated... maybe it's prime? Do a fast check.
|
||||
if (baillie_psw_probable_prime(n)) {
|
||||
return n;
|
||||
}
|
||||
|
||||
[[nodiscard]] static consteval bool is_prime(std::size_t n) { return (n > 1) && find_first_factor(n) == n; }
|
||||
};
|
||||
// If we're here, we know `n` is composite, so continue with trial division for all odd numbers.
|
||||
std::uintmax_t factor = first_100_primes.back() + 2u;
|
||||
while (factor * factor <= n) {
|
||||
if (n % factor == 0u) {
|
||||
return factor;
|
||||
}
|
||||
factor += 2u;
|
||||
}
|
||||
|
||||
return n; // Technically unreachable.
|
||||
}
|
||||
|
||||
} // namespace mp_units::detail
|
||||
|
||||
@ -51,13 +51,6 @@ import std;
|
||||
|
||||
namespace mp_units {
|
||||
|
||||
namespace detail {
|
||||
|
||||
// Higher numbers use fewer trial divisions, at the price of more storage space.
|
||||
using factorizer = wheel_factorizer<4>;
|
||||
|
||||
} // namespace detail
|
||||
|
||||
#if defined MP_UNITS_COMP_CLANG || MP_UNITS_COMP_CLANG < 18
|
||||
|
||||
MP_UNITS_EXPORT template<symbol_text Symbol>
|
||||
@ -629,14 +622,10 @@ using common_magnitude_type = decltype(_common_magnitude_type_impl(M));
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// `mag()` implementation.
|
||||
|
||||
// Sometimes we need to give the compiler a "shortcut" when factorizing large numbers (specifically, numbers whose
|
||||
// _first factor_ is very large). If we don't, we can run into limits on the number of constexpr steps or iterations.
|
||||
//
|
||||
// To provide the first factor for a given number, specialize this variable template.
|
||||
//
|
||||
// WARNING: The program behaviour will be undefined if you provide a wrong answer, so check your math!
|
||||
MP_UNITS_EXPORT template<std::intmax_t N>
|
||||
constexpr std::optional<std::intmax_t> known_first_factor = std::nullopt;
|
||||
[[deprecated("`known_first_factor` is no longer necessary and can simply be removed")]]
|
||||
constexpr std::optional<std::intmax_t>
|
||||
known_first_factor = std::nullopt;
|
||||
|
||||
namespace detail {
|
||||
|
||||
@ -646,12 +635,7 @@ template<std::intmax_t N>
|
||||
struct prime_factorization {
|
||||
[[nodiscard]] static consteval std::intmax_t get_or_compute_first_factor()
|
||||
{
|
||||
constexpr auto opt = known_first_factor<N>;
|
||||
if constexpr (opt.has_value()) {
|
||||
return opt.value(); // NOLINT(bugprone-unchecked-optional-access)
|
||||
} else {
|
||||
return static_cast<std::intmax_t>(factorizer::find_first_factor(N));
|
||||
}
|
||||
return static_cast<std::intmax_t>(find_first_factor(N));
|
||||
}
|
||||
|
||||
static constexpr std::intmax_t first_base = get_or_compute_first_factor();
|
||||
|
||||
@ -35,9 +35,6 @@
|
||||
MP_UNITS_EXPORT
|
||||
namespace mp_units {
|
||||
|
||||
template<>
|
||||
constexpr std::optional<std::intmax_t> known_first_factor<334'524'384'739> = 334'524'384'739;
|
||||
|
||||
namespace hep {
|
||||
|
||||
// energy
|
||||
|
||||
@ -31,9 +31,6 @@ import mp_units;
|
||||
using namespace units;
|
||||
using namespace units::detail;
|
||||
|
||||
template<>
|
||||
constexpr std::optional<std::intmax_t> units::known_first_factor<9223372036854775783> = 9223372036854775783;
|
||||
|
||||
namespace {
|
||||
|
||||
// A set of non-standard bases for testing purposes.
|
||||
@ -182,20 +179,6 @@ static_assert(std::is_same_v<decltype(get_base(power_v<mag_2, 5, 8>{})), mag_2_>
|
||||
// mag_ratio<16'605'390'666'050, 10'000'000'000'000>();
|
||||
// }
|
||||
|
||||
// SECTION("Can bypass computing primes by providing known_first_factor<N>")
|
||||
// {
|
||||
// // Sometimes, even wheel factorization isn't enough to handle the compilers' limits on constexpr steps and/or
|
||||
// // iterations. To work around these cases, we can explicitly provide the correct answer directly to the
|
||||
// compiler.
|
||||
// //
|
||||
// // In this case, we test that we can represent the largest prime that fits in a signed 64-bit int. The reason
|
||||
// this
|
||||
// // test can pass is that we have provided the answer, by specializing the `known_first_factor` variable template
|
||||
// // above in this file.
|
||||
// mag<9'223'372'036'854'775'783>();
|
||||
// }
|
||||
// }
|
||||
|
||||
// TEST_CASE("magnitude converts to numerical value")
|
||||
// {
|
||||
// SECTION("Positive integer powers of integer bases give integer values")
|
||||
|
||||
@ -35,51 +35,6 @@ namespace {
|
||||
|
||||
inline constexpr auto MAX_U64 = std::numeric_limits<std::uint64_t>::max();
|
||||
|
||||
template<std::size_t BasisSize, std::size_t... Is>
|
||||
constexpr bool check_primes(std::index_sequence<Is...>)
|
||||
{
|
||||
return ((Is < 2 || wheel_factorizer<BasisSize>::is_prime(Is) == is_prime_by_trial_division(Is)) && ...);
|
||||
}
|
||||
|
||||
static_assert(check_primes<2>(std::make_index_sequence<122>{}));
|
||||
|
||||
// This is the smallest number that can catch the bug where we use only _prime_ numbers in the first wheel, rather than
|
||||
// numbers which are _coprime to the basis_.
|
||||
//
|
||||
// The basis for N = 4 is {2, 3, 5, 7}, so the wheel size is 210. 11 * 11 = 121 is within the first wheel. It is
|
||||
// coprime with every element of the basis, but it is _not_ prime. If we keep only prime numbers, then we will neglect
|
||||
// using numbers of the form (210 * n + 121) as trial divisors, which is a problem if any are prime. For n = 1, we have
|
||||
// a divisor of (210 + 121 = 331), which happens to be prime but will not be used. Thus, (331 * 331 = 109561) is a
|
||||
// composite number which could wrongly appear prime if we skip over 331.
|
||||
static_assert(wheel_factorizer<4>::is_prime(109'561) == is_prime_by_trial_division(109'561));
|
||||
|
||||
static_assert(wheel_factorizer<1>::coprimes_in_first_wheel.size() == 1);
|
||||
static_assert(wheel_factorizer<2>::coprimes_in_first_wheel.size() == 2);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel.size() == 8);
|
||||
static_assert(wheel_factorizer<4>::coprimes_in_first_wheel.size() == 48);
|
||||
static_assert(wheel_factorizer<5>::coprimes_in_first_wheel.size() == 480);
|
||||
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[0] == 1);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[1] == 7);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[2] == 11);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[3] == 13);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[4] == 17);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[5] == 19);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[6] == 23);
|
||||
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[7] == 29);
|
||||
|
||||
static_assert(!wheel_factorizer<1>::is_prime(0));
|
||||
static_assert(!wheel_factorizer<1>::is_prime(1));
|
||||
static_assert(wheel_factorizer<1>::is_prime(2));
|
||||
|
||||
static_assert(!wheel_factorizer<2>::is_prime(0));
|
||||
static_assert(!wheel_factorizer<2>::is_prime(1));
|
||||
static_assert(wheel_factorizer<2>::is_prime(2));
|
||||
|
||||
static_assert(!wheel_factorizer<3>::is_prime(0));
|
||||
static_assert(!wheel_factorizer<3>::is_prime(1));
|
||||
static_assert(wheel_factorizer<3>::is_prime(2));
|
||||
|
||||
// Modular arithmetic.
|
||||
static_assert(add_mod(1u, 2u, 5u) == 3u);
|
||||
static_assert(add_mod(4u, 4u, 5u) == 3u);
|
||||
@ -122,4 +77,103 @@ static_assert(miller_rabin_probable_prime(2u, 9'007'199'254'740'881u), "Large kn
|
||||
|
||||
static_assert(miller_rabin_probable_prime(2u, 18'446'744'073'709'551'557u), "Largest 64-bit prime");
|
||||
|
||||
// Jacobi symbols --- a building block for the Strong Lucas probable prime test, needed for Baillie-PSW.
|
||||
static_assert(jacobi_symbol(1, 1u) == 1, "Jacobi symbol always 1 when 'numerator' is 1");
|
||||
static_assert(jacobi_symbol(1, 3u) == 1, "Jacobi symbol always 1 when 'numerator' is 1");
|
||||
static_assert(jacobi_symbol(1, 5u) == 1, "Jacobi symbol always 1 when 'numerator' is 1");
|
||||
static_assert(jacobi_symbol(1, 987654321u) == 1, "Jacobi symbol always 1 when 'numerator' is 1");
|
||||
|
||||
static_assert(jacobi_symbol(3, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1");
|
||||
static_assert(jacobi_symbol(5, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1");
|
||||
static_assert(jacobi_symbol(-1234567890, 1u) == 1, "Jacobi symbol always 1 when 'denominator' is 1");
|
||||
|
||||
static_assert(jacobi_symbol(10, 5u) == 0, "Jacobi symbol always 0 when there's a common factor");
|
||||
static_assert(jacobi_symbol(25, 15u) == 0, "Jacobi symbol always 0 when there's a common factor");
|
||||
static_assert(jacobi_symbol(-24, 9u) == 0, "Jacobi symbol always 0 when there's a common factor");
|
||||
|
||||
static_assert(jacobi_symbol(14, 9u) == +jacobi_symbol(7, 9u),
|
||||
"Divide numerator by 2: positive when (denom % 8) in {1, 7}");
|
||||
static_assert(jacobi_symbol(14, 15u) == +jacobi_symbol(7, 15u),
|
||||
"Divide numerator by 2: positive when (denom % 8) in {1, 7}");
|
||||
static_assert(jacobi_symbol(14, 11u) == -jacobi_symbol(7, 11u),
|
||||
"Divide numerator by 2: negative when (denom % 8) in {3, 5}");
|
||||
static_assert(jacobi_symbol(14, 13u) == -jacobi_symbol(7, 13u),
|
||||
"Divide numerator by 2: negative when (denom % 8) in {3, 5}");
|
||||
|
||||
static_assert(jacobi_symbol(19, 9u) == +jacobi_symbol(9, 19u), "Flip is identity when (n % 4) = 1");
|
||||
static_assert(jacobi_symbol(17, 7u) == +jacobi_symbol(7, 17u), "Flip is identity when (a % 4) = 1");
|
||||
static_assert(jacobi_symbol(19, 7u) == -jacobi_symbol(9, 7u), "Flip changes sign when (n % 4) = 3 and (a % 4) = 3");
|
||||
|
||||
static_assert(jacobi_symbol(1001, 9907u) == -1, "Example from Wikipedia page");
|
||||
static_assert(jacobi_symbol(19, 45u) == 1, "Example from Wikipedia page");
|
||||
static_assert(jacobi_symbol(8, 21u) == -1, "Example from Wikipedia page");
|
||||
static_assert(jacobi_symbol(5, 21u) == 1, "Example from Wikipedia page");
|
||||
|
||||
// Tests for perfect square finder
|
||||
static_assert(is_perfect_square(0u));
|
||||
static_assert(is_perfect_square(1u));
|
||||
static_assert(!is_perfect_square(2u));
|
||||
static_assert(is_perfect_square(4u));
|
||||
|
||||
constexpr std::uint64_t BIG_SQUARE = [](auto x) { return x * x; }((std::uint64_t{1u} << 32) - 1u);
|
||||
static_assert(!is_perfect_square(BIG_SQUARE - 1u));
|
||||
static_assert(is_perfect_square(BIG_SQUARE));
|
||||
static_assert(!is_perfect_square(BIG_SQUARE + 1u));
|
||||
|
||||
// Tests for the Strong Lucas Probable Prime test.
|
||||
static_assert(as_int(LucasDParameter{.mag = 5, .pos = true}) == 5);
|
||||
static_assert(as_int(LucasDParameter{.mag = 7, .pos = false}) == -7);
|
||||
|
||||
static_assert(as_int(LucasDParameter{}) == 5, "First D parameter in the sequence is 5");
|
||||
static_assert(as_int(successor(LucasDParameter{})) == -7, "Incrementing adds 2 to the mag, and flips the sign");
|
||||
static_assert(as_int(successor(successor(LucasDParameter{}))) == 9);
|
||||
static_assert(as_int(successor(successor(successor(LucasDParameter{})))) == -11);
|
||||
|
||||
static_assert(strong_lucas_probable_prime(3u), "Known small prime");
|
||||
static_assert(strong_lucas_probable_prime(5u), "Known small prime");
|
||||
static_assert(strong_lucas_probable_prime(7u), "Known small prime");
|
||||
static_assert(!strong_lucas_probable_prime(9u), "Known small composite");
|
||||
|
||||
// Test some Miller-Rabin pseudoprimes (https://oeis.org/A001262), which should NOT be marked prime.
|
||||
static_assert(!strong_lucas_probable_prime(2047u), "Miller-Rabin pseudoprime");
|
||||
static_assert(!strong_lucas_probable_prime(3277u), "Miller-Rabin pseudoprime");
|
||||
static_assert(!strong_lucas_probable_prime(486737u), "Miller-Rabin pseudoprime");
|
||||
|
||||
// Test some Strong Lucas pseudoprimes (https://oeis.org/A217255).
|
||||
static_assert(strong_lucas_probable_prime(5459u), "Strong Lucas pseudoprime");
|
||||
static_assert(strong_lucas_probable_prime(5777u), "Strong Lucas pseudoprime");
|
||||
static_assert(strong_lucas_probable_prime(10877u), "Strong Lucas pseudoprime");
|
||||
static_assert(strong_lucas_probable_prime(324899u), "Strong Lucas pseudoprime");
|
||||
|
||||
// Test some actual primes
|
||||
static_assert(strong_lucas_probable_prime(225'653'407'801u), "Large known prime");
|
||||
static_assert(strong_lucas_probable_prime(334'524'384'739u), "Large known prime");
|
||||
static_assert(strong_lucas_probable_prime(9'007'199'254'740'881u), "Large known prime");
|
||||
|
||||
static_assert(strong_lucas_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime");
|
||||
|
||||
// Tests for Baillie-PSW, which is known to be correct for all 64-bit integers.
|
||||
static_assert(baillie_psw_probable_prime(3u), "Known small prime");
|
||||
static_assert(baillie_psw_probable_prime(5u), "Known small prime");
|
||||
static_assert(baillie_psw_probable_prime(7u), "Known small prime");
|
||||
static_assert(!baillie_psw_probable_prime(9u), "Known small composite");
|
||||
|
||||
// Test some Miller-Rabin pseudoprimes (https://oeis.org/A001262), which should NOT be marked prime.
|
||||
static_assert(!baillie_psw_probable_prime(2047u), "Miller-Rabin pseudoprime");
|
||||
static_assert(!baillie_psw_probable_prime(3277u), "Miller-Rabin pseudoprime");
|
||||
static_assert(!baillie_psw_probable_prime(486737u), "Miller-Rabin pseudoprime");
|
||||
|
||||
// Test some Strong Lucas pseudoprimes (https://oeis.org/A217255), which should NOT be marked prime.
|
||||
static_assert(!baillie_psw_probable_prime(5459u), "Strong Lucas pseudoprime");
|
||||
static_assert(!baillie_psw_probable_prime(5777u), "Strong Lucas pseudoprime");
|
||||
static_assert(!baillie_psw_probable_prime(10877u), "Strong Lucas pseudoprime");
|
||||
static_assert(!baillie_psw_probable_prime(324899u), "Strong Lucas pseudoprime");
|
||||
|
||||
// Test some actual primes
|
||||
static_assert(baillie_psw_probable_prime(225'653'407'801u), "Large known prime");
|
||||
static_assert(baillie_psw_probable_prime(334'524'384'739u), "Large known prime");
|
||||
static_assert(baillie_psw_probable_prime(9'007'199'254'740'881u), "Large known prime");
|
||||
|
||||
static_assert(baillie_psw_probable_prime(18'446'744'073'709'551'557u), "Largest 64-bit prime");
|
||||
|
||||
} // namespace
|
||||
|
||||
Reference in New Issue
Block a user