Skip to content

Commit

Permalink
update math_functions.hpp and is_prime.hpp
Browse files Browse the repository at this point in the history
  • Loading branch information
@i80287
i80287 committed Nov 12, 2024
1 parent d20e812 commit 9295bd5
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 39 deletions.
24 changes: 12 additions & 12 deletions number_theory/is_prime.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,9 +91,9 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_prp(uint64_t n, uint64_t a) noexce
* the Jacobi symbol]
**********************************************************************************************/
template <bool BasicChecks = true>
ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_lucas_prp(uint64_t n, uint32_t p,
ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_lucas_prp(uint64_t n, uint16_t p,
int32_t q) noexcept {
const int64_t d = static_cast<int64_t>(uint64_t{p} * uint64_t{p}) - int64_t{q} * 4;
const int64_t d = int64_t{uint32_t{p} * uint32_t{p}} - int64_t{q} * 4;
if constexpr (BasicChecks) {
/* Check if p*p - 4*q == 0. */
if (unlikely(d == 0)) {
Expand Down Expand Up @@ -139,16 +139,15 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_lucas_prp(uint64_t n, uint32_t p,
* make sure U_s == 0 mod n or V_((2^t)*s) == 0 mod n,
* for some t, 0 <= t < r
*/
uint64_t uh = 1; // Initial value for U_1
uint64_t vl = 2; // Initial value for V_0
uint64_t vh = p; // Initial value for V_1
uint64_t uh = 1; // Initial value for U_1
uint64_t vl = 2; // Initial value for V_0
uint64_t vh = uint64_t{p}; // Initial value for V_1
uint64_t ql = 1;
uint64_t qh = 1;
// q mod n
const uint64_t widen_q =
(q >= 0 ? static_cast<uint32_t>(q)
: (n - (static_cast<uint64_t>(-static_cast<uint32_t>(q)) % n))) %
n;
const uint64_t widen_q = (q >= 0 ? static_cast<uint32_t>(q)
: (n - static_cast<uint64_t>(-static_cast<uint32_t>(q)) % n)) %
n;
CONFIG_ASSUME_STATEMENT(widen_q < n);
// n >= 3 => n - 1 >= 2 => n - 1 >= 1 => s >= 1
for (uint32_t j = ::math_functions::log2_floor(s); j != 0; j--) {
Expand Down Expand Up @@ -290,7 +289,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_lucas_prp(uint64_t n, uint32_t p,
tmp_vl = vl_vl - ql_2;
vl = vl_vl >= ql_2 ? tmp_vl : tmp_vl + n;
CONFIG_ASSUME_STATEMENT(vl < n);
CONFIG_ASSUME_STATEMENT(vl == (uint128_t(n) + vl_vl - ql_2) % n);
CONFIG_ASSUME_STATEMENT(vl == (uint128_t{n} + vl_vl - ql_2) % n);

if (vl == 0) {
return true;
Expand Down Expand Up @@ -355,6 +354,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_selfridge_prp(uint64_t n) noexcept
break;
case -1: {
CONFIG_ASSUME_STATEMENT(d <= kMaxD + kStep * 2);
CONFIG_ASSUME_STATEMENT(-kMaxD - kStep <= d);
CONFIG_ASSUME_STATEMENT((1 - d) % 4 == 0);
const std::int32_t q = (1 - d) / 4;
CONFIG_ASSUME_STATEMENT(1 - 4 * q == d);
Expand Down Expand Up @@ -586,8 +586,8 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_selfridge_prp(uint64_t n) noexcept
}

[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR bool is_mersenne_prime(const uint128_t n) noexcept {
const auto np1 = n + 1;
if (!is_power_of_two(np1)) {
const uint128_t np1 = n + 1;
if (!::math_functions::is_power_of_two(np1)) {
return false;
}
const auto [q, p] = ::math_functions::extract_pow2(np1);
Expand Down
72 changes: 45 additions & 27 deletions number_theory/math_functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR uint64_t bin_pow_mod(uint64_t n, uint64_t p, uint
return y;
}

[[nodiscard]] ATTRIBUTE_CONST constexpr uint32_t isqrt(uint64_t n) noexcept {
[[nodiscard]] ATTRIBUTE_CONST constexpr uint32_t isqrt(const uint64_t n) noexcept {
/**
* In the runtime `sqrtl` is used (but not for the msvc prior to the c++20).
*/
Expand Down Expand Up @@ -236,7 +236,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR uint64_t bin_pow_mod(uint64_t n, uint64_t p, uint

#if defined(INTEGERS_128_BIT_HPP)

[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint64_t isqrt(uint128_t n) noexcept {
[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint64_t isqrt(const uint128_t n) noexcept {
/**
* See Hackers Delight Chapter 11.
*/
Expand Down Expand Up @@ -282,7 +282,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR uint64_t bin_pow_mod(uint64_t n, uint64_t p, uint
* `uint32_t(std::cbrt(3375.0))` may be equal to 14
*/

#if defined(__GNUG__) && !defined(__clang__)
#if defined(__GNUG__) && !defined(__clang__) && CONFIG_HAS_AT_LEAST_CXX_17
[[maybe_unused]] const auto n_original_value = n;
#endif

Expand Down Expand Up @@ -341,7 +341,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR uint64_t bin_pow_mod(uint64_t n, uint64_t p, uint
/// @note ⌊n^0.25⌋ = ⌊⌊n^0.5⌋^0.5⌋ (see Hackers Delight Chapter 11, ex.1)
/// @param[in] n
/// @return ⌊n^0.25⌋
[[nodiscard]] ATTRIBUTE_CONST constexpr uint32_t ifrrt(uint64_t n) noexcept {
[[nodiscard]] ATTRIBUTE_CONST constexpr uint32_t ifrrt(const uint64_t n) noexcept {
return ::math_functions::isqrt(::math_functions::isqrt(n));
}

Expand All @@ -351,7 +351,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR uint64_t bin_pow_mod(uint64_t n, uint64_t p, uint
/// It can be shown that ⌊n^0.25⌋ = ⌊⌊n^0.5⌋^0.5⌋
/// @param[in] n
/// @return
[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint32_t ifrrt(uint128_t n) noexcept {
[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint32_t ifrrt(const uint128_t n) noexcept {
return ::math_functions::isqrt(::math_functions::isqrt(n));
}

Expand Down Expand Up @@ -468,7 +468,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR
/// @brief This function reverses bits of the @a `b`
/// @param[in] b
/// @return 8-bit number whose bits are reversed bits of the @a `b`.
[[nodiscard]] ATTRIBUTE_CONST constexpr uint8_t bit_reverse(uint8_t b) noexcept {
[[nodiscard]] ATTRIBUTE_CONST constexpr uint8_t bit_reverse(const uint8_t b) noexcept {
// See https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2
return static_cast<uint8_t>(((b * 0x80200802ULL) & 0x0884422110ULL) * 0x0101010101ULL >> 32U);
}
Expand Down Expand Up @@ -1413,7 +1413,9 @@ ATTRIBUTE_ALWAYS_INLINE ATTRIBUTE_CONST constexpr uint32_t base_b_len_impl(
template <typename T>
[[nodiscard]] ATTRIBUTE_CONST constexpr uint32_t base_b_len(T value,
const uint8_t base = 10) noexcept {
static_assert(::math_functions::detail::is_integral_v<T>);
static_assert(::math_functions::detail::is_integral_v<T> && !std::is_same_v<T, bool> &&
!std::is_same_v<T, char>,
"integral type (not bool or char) expected in the base_b_len");

if constexpr (::math_functions::detail::is_signed_v<T>) {
const uint32_t is_negative = uint32_t{value < 0};
Expand All @@ -1424,7 +1426,7 @@ template <typename T>
}
}

/// @brief For n > 0 returns ⌈log_2(n)⌉. For n = 0 returns (uint32_t)-1
/// @brief For n > 0 returns ⌊log_2(n)⌋. For n = 0 returns (uint32_t)-1
/// @tparam UIntType unsigned integral type (at least unsigned int in size)
/// @param[in] n
/// @return
Expand Down Expand Up @@ -1461,7 +1463,7 @@ template <class UIntType>
#endif
[[nodiscard]]
ATTRIBUTE_ALWAYS_INLINE ATTRIBUTE_CONST constexpr uint32_t log2_ceil(const UIntType n) noexcept {
return ::math_functions::log2_floor(n) + ((n & (n - 1)) != 0);
return ::math_functions::log2_floor(n) + uint32_t{(n & (n - 1)) != 0};
}

template <class UIntType>
Expand Down Expand Up @@ -1680,29 +1682,15 @@ template <typename T>
return (x || y) && (y || z) && (x || z);
}

template <class T>
#if CONFIG_HAS_CONCEPTS

template <::math_functions::detail::unsigned_integral T>
requires ::math_functions::detail::unsigned_integral<T>
#endif
[[nodiscard]] ATTRIBUTE_CONST constexpr T next_even(T n) noexcept {
static_assert(::math_functions::detail::is_unsigned_v<T>, "unsigned integral type expected");
return n + 2 - n % 2;
}

#else

// clang-format off

template <class T>
[[nodiscard]]
ATTRIBUTE_CONST
constexpr
std::enable_if_t<::math_functions::detail::is_unsigned_v<T>, T> next_even(T n) noexcept {
return n + 2 - n % 2;
}

// clang-format on

#endif

template <class FloatType>
struct SumSinCos {
FloatType sines_sum;
Expand Down Expand Up @@ -2923,6 +2911,36 @@ template <class T>
return ::math_functions::arange(T{0}, n);
}

/// @brief Return vector of elements {log2(0), log2(1), log2(2), log2(3), ..., log2(n)}
/// @note Here log2(0) := -1
/// @param n
/// @return
[[nodiscard]] CONSTEXPR_VECTOR std::vector<uint32_t> log2_arange(const uint32_t n) {
std::vector<uint32_t> values(size_t{n} + 1);
values[0] = static_cast<uint32_t>(-1);
for (size_t i = 1; i <= n; i++) {
values[i] = values[i / 2] + 1;
}

return values;
}

/// @brief Return vector of elements {0! mod m, 1! mod m, 2! mod m, 3! mod m, ..., n! mod m}
/// @param n
/// @return
[[nodiscard]]
CONSTEXPR_VECTOR std::vector<uint32_t> factorial_arange_mod_m(const uint32_t n, const uint32_t m) {
std::vector<uint32_t> values(size_t{n} + 1);
uint32_t current_factorial = m != 1 ? 1u : 0u;
values[0] = current_factorial;
for (size_t i = 1; i <= n; i++) {
current_factorial = static_cast<uint32_t>((uint64_t{current_factorial} * uint64_t{i}) % m);
values[i] = current_factorial;
}

return values;
}

} // namespace math_functions

// NOLINTEND(cppcoreguidelines-avoid-magic-numbers)
Expand Down
32 changes: 32 additions & 0 deletions number_theory/test_math_functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3471,6 +3471,36 @@ void test_arange() {
assert((arange(11, 11, 11).empty()));
}

void test_log2_arange() {
log_tests_started();

using std::vector;

assert(log2_arange(0) == vector<uint32_t>{static_cast<uint32_t>(-1)});
const uint32_t n = 1'000'000;
const vector<uint32_t> range = log2_arange(n);
assert(range.size() == n + 1);
for (uint32_t i = 0; i <= n; i++) {
assert(range[i] == log2_floor(i));
}
}

void test_factorial_arange_mod_m() {
log_tests_started();

for (const uint32_t m : {2u, 4u, static_cast<uint32_t>(1e7) + 9}) {
for (const uint32_t n : {10u, 1000u, 100000u}) {
const std::vector<uint32_t> fact_range = factorial_arange_mod_m(n, m);
uint32_t factorial = 1;
assert(fact_range[0] == factorial);
for (uint32_t i = 1; size_t{i} <= size_t{n}; i++) {
factorial = static_cast<uint32_t>((uint64_t{factorial} * uint64_t{i}) % m);
assert(fact_range[i] == factorial);
}
}
}
}

void test_masked_popcount_sum() noexcept {
log_tests_started();

Expand Down Expand Up @@ -3523,5 +3553,7 @@ int main() {
test_solve_factorial_congruence();
test_powers_sum();
test_arange();
test_log2_arange();
test_factorial_arange_mod_m();
test_masked_popcount_sum();
}

0 comments on commit 9295bd5

Please sign in to comment.