Skip to content

Commit

Permalink
update gcd() in the math_functions.hpp
Browse files Browse the repository at this point in the history
  • Loading branch information
@i80287
i80287 committed Nov 13, 2024
1 parent 20d61b5 commit 67a3cda
Show file tree
Hide file tree
Showing 3 changed files with 191 additions and 74 deletions.
2 changes: 1 addition & 1 deletion number_theory/is_prime.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ ATTRIBUTE_CONST I128_CONSTEXPR bool is_strong_lucas_prp(uint64_t n, uint16_t p,

// NOLINTNEXTLINE(bugprone-implicit-widening-of-multiplication-result)
const int128_t rhs = int128_t{int64_t{2} * int64_t{q}} * int128_t{d};
if (unlikely(std::gcd(n, rhs) != 1)) {
if (unlikely(math_functions::gcd(n, rhs) != 1)) {
// is_strong_lucas_prp requires gcd(n, 2 * q * (p * p - 4 * q)) = 1
return false;
}
Expand Down
146 changes: 114 additions & 32 deletions number_theory/math_functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3043,22 +3043,15 @@ ATTRIBUTE_NODISCARD CONSTEXPR_VECTOR auto weighted_median(Range&& range ATTRIBUT

#endif

} // namespace math_functions

// NOLINTEND(cppcoreguidelines-avoid-magic-numbers)

#if defined(INTEGERS_128_BIT_HPP)

namespace std {
#ifdef INTEGERS_128_BIT_HPP

// NOLINTBEGIN(cert-dcl58-cpp)
namespace detail {

/// @brief Computes greaters common divisor of @a `a` and @a `b`
/// using Stein's algorithm (binary gcd). Here gcd(0, 0) = 0.
/// @param[in] a
/// @param[in] b
/// @return gcd(a, b)
[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint128_t gcd(uint128_t a, uint128_t b) noexcept {
ATTRIBUTE_NODISCARD
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(uint128_t a, uint128_t b) noexcept {
if (unlikely(a == 0)) {
return b;
}
Expand Down Expand Up @@ -3088,37 +3081,126 @@ namespace std {
}
}

[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint128_t gcd(uint64_t a, int128_t b) noexcept {
const uint128_t b0 = ::math_functions::uabs(b);
if (unlikely(b0 == 0)) {
return a;
ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(uint128_t a, int128_t b) noexcept {
return ::math_functions::detail::gcd(a, ::math_functions::uabs(b));
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(int128_t a, uint128_t b) noexcept {
return ::math_functions::detail::gcd(::math_functions::uabs(a), b);
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR int128_t gcd(int128_t a, int128_t b) noexcept {
const int128_t value =
static_cast<int128_t>(::math_functions::detail::gcd(::math_functions::uabs(a), b));
CONFIG_ASSUME_STATEMENT(value >= 0);
return value;
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(uint128_t a, uint64_t b) noexcept {
if ((config::is_constant_evaluated() || config::is_gcc_constant_p(a <= b)) && a <= b) {
return std::gcd(static_cast<uint64_t>(a), b);
}

// gcd(a, b) = gcd(a, b0) = gcd(b0, a % b0) = gcd(a1, b1)
const uint128_t a1 = b0;
// b1 = a % b0
const uint64_t b1 = a < b0 ? a : a % static_cast<uint64_t>(b0); // a < 2^64 => b1 < 2^64
if (b1 == 0) {
return a1;
if (unlikely(b == 0)) {
return a;
}
// gcd(a1, b1) = gcd(b1, a1 % b1) = gcd(a2, b2)
const uint64_t a2 = b1; // b1 < 2^64 => a2 < 2^64
// b2 = a1 % b1
// a1 = b0, b1 = a % b0 => b1 < a1
const uint64_t b2 = static_cast<uint64_t>(a1 % b1); // b1 < 2^64 => b2 = a1 % b1 < 2^64
return std::gcd(a2, b2);
// gcd(a, b) = gcd(b, a % b) = gcd(a % b, b)
return std::gcd(static_cast<uint64_t>(a % b), b);
}

[[nodiscard]] ATTRIBUTE_CONST I128_CONSTEXPR uint128_t gcd(int128_t a, uint64_t b) noexcept {
return std::gcd(b, a);
ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(uint64_t a, uint128_t b) noexcept {
return ::math_functions::detail::gcd(b, a);
}

// NOLINTEND(cert-dcl58-cpp)
ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(uint128_t a, int64_t b) noexcept {
return ::math_functions::detail::gcd(a, ::math_functions::uabs(b));
}

} // namespace std
ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR uint128_t gcd(int64_t a, uint128_t b) noexcept {
return ::math_functions::detail::gcd(b, a);
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR int128_t gcd(uint64_t a, int128_t b) noexcept {
const int128_t value =
static_cast<int128_t>(::math_functions::detail::gcd(a, ::math_functions::uabs(b)));
CONFIG_ASSUME_STATEMENT(value >= 0);
return value;
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR int128_t gcd(int128_t a, uint64_t b) noexcept {
return ::math_functions::detail::gcd(b, a);
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR int128_t gcd(int128_t a, int64_t b) noexcept {
return ::math_functions::detail::gcd(a, ::math_functions::uabs(b));
}

ATTRIBUTE_NODISCARD
ATTRIBUTE_ALWAYS_INLINE
ATTRIBUTE_CONST
I128_CONSTEXPR int128_t gcd(int64_t a, int128_t b) noexcept {
return ::math_functions::detail::gcd(b, a);
}

} // namespace detail

#endif // INTEGERS_128_BIT_HPP

/// @brief Computes greaters common divisor of @a `a` and @a `b`
/// using Stein's algorithm (binary gcd). Here gcd(0, 0) = 0.
/// @param[in] a
/// @param[in] b
/// @return gcd(a, b)
template <class M, class N>
[[nodiscard]]
ATTRIBUTE_ALWAYS_INLINE ATTRIBUTE_CONST constexpr std::common_type_t<M, N> gcd(M m, N n) noexcept {
static_assert(
::math_functions::detail::is_integral_v<M> && ::math_functions::detail::is_integral_v<N>,
"math_functions::gcd arguments must be integers");

#if defined(INTEGERS_128_BIT_HPP)
if constexpr (sizeof(M) <= sizeof(uint64_t) && sizeof(N) <= sizeof(uint64_t)) {
#endif
return std::gcd(m, n);
#if defined(INTEGERS_128_BIT_HPP)
} else {
return ::math_functions::detail::gcd(m, n);
}
#endif
}

} // namespace math_functions

#undef CONSTEXPR_VECTOR
#ifdef MATH_FUNCTIONS_HAS_BIT
#undef MATH_FUNCTIONS_HAS_BIT
Expand Down
117 changes: 76 additions & 41 deletions number_theory/test_math_functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@
using namespace math_functions;
// NOLINTNEXTLINE(google-build-using-namespace)
using namespace test_tools;
using std::gcd;
using std::size_t;
using std::uint32_t;
using std::uint64_t;
Expand Down Expand Up @@ -1869,52 +1868,88 @@ void test_general_asserts() noexcept {
I128_ASSERT_THAT(base_b_len(-int128_t{101}) == 4);
I128_ASSERT_THAT(base_b_len(static_cast<int128_t>(uint128_t{1} << 127U)) == 40);

I128_ASSERT_THAT(gcd(uint128_t{1}, uint128_t{1}) == 1);
I128_ASSERT_THAT(gcd(uint128_t{3}, uint128_t{7}) == 1);
I128_ASSERT_THAT(gcd(uint128_t{0}, uint128_t{112378432}) == 112378432);
I128_ASSERT_THAT(gcd(uint128_t{112378432}, uint128_t{0}) == 112378432);
I128_ASSERT_THAT(gcd(uint128_t{429384832}, uint128_t{324884}) == 4);
I128_ASSERT_THAT(gcd(uint128_t{18446744073709551521ULL}, uint128_t{18446744073709551533ULL}) ==
1);
I128_ASSERT_THAT(gcd(uint128_t{18446744073709551521ULL} * 18446744073709551521ULL,
uint128_t{18446744073709551521ULL}) == 18446744073709551521ULL);
I128_ASSERT_THAT(gcd(uint128_t{23999993441ULL} * 23999993377ULL,
uint128_t{23999992931ULL} * 23999539633ULL) == 1);
I128_ASSERT_THAT(gcd(uint128_t{2146514599U} * 2146514603U * 2146514611U,
uint128_t{2146514611U} * 2146514621U * 2146514647U) == 2146514611ULL);
I128_ASSERT_THAT(gcd(uint128_t{2146514599U} * 2146514603U * 2146514611U * 2,
uint128_t{2146514599U} * 2146514603U * 2146514611U * 3) ==
uint128_t{2146514599U} * 2146514603U * 2146514611U);
I128_ASSERT_THAT(gcd(uint128_t{100000000000000003ULL} * 1000000000000000003ULL,
uint128_t{1000000000000000003ULL} * 1000000000000000009ULL) ==
1000000000000000003ULL);
I128_ASSERT_THAT(gcd(uint128_t{3} * 2 * 5 * 7 * 11 * 13 * 17 * 19,
uint128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(gcd(uint128_t{1000000000000000009ULL},
uint128_t{1000000000000000009ULL} * 1000000000000000009ULL) ==
1000000000000000009ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{1}, uint128_t{1}) == 1);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{3}, uint128_t{7}) == 1);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{0}, uint128_t{112378432}) == 112378432);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{112378432}, uint128_t{0}) == 112378432);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{429384832}, uint128_t{324884}) == 4);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551521ULL},
uint128_t{18446744073709551533ULL}) == 1);
I128_ASSERT_THAT(
gcd(uint128_t{0}, uint128_t{1000000000000000009ULL} * 1000000000000000009ULL) ==
uint128_t{1000000000000000009ULL} * 1000000000000000009ULL);
I128_ASSERT_THAT(gcd(uint128_t{18446744073709551557ULL}, uint128_t{0}) ==
math_functions::gcd(uint128_t{18446744073709551521ULL} * 18446744073709551521ULL,
uint128_t{18446744073709551521ULL}) == 18446744073709551521ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{23999993441ULL} * 23999993377ULL,
uint128_t{23999992931ULL} * 23999539633ULL) == 1);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{2146514599U} * 2146514603U * 2146514611U,
uint128_t{2146514611U} * 2146514621U * 2146514647U) ==
2146514611ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{2146514599U} * 2146514603U * 2146514611U * 2,
uint128_t{2146514599U} * 2146514603U * 2146514611U * 3) ==
uint128_t{2146514599U} * 2146514603U * 2146514611U);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{100000000000000003ULL} * 1000000000000000003ULL,
uint128_t{1000000000000000003ULL} *
1000000000000000009ULL) == 1000000000000000003ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{3} * 2 * 5 * 7 * 11 * 13 * 17 * 19,
uint128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{1000000000000000009ULL},
uint128_t{1000000000000000009ULL} *
1000000000000000009ULL) == 1000000000000000009ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{0}, uint128_t{1000000000000000009ULL} *
1000000000000000009ULL) ==
uint128_t{1000000000000000009ULL} * 1000000000000000009ULL);
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551557ULL}, uint128_t{0}) ==
18446744073709551557ULL);

I128_ASSERT_THAT(gcd(uint64_t{2}, int128_t{4}) == 2);
I128_ASSERT_THAT(gcd(uint64_t{2}, int128_t{-4}) == 2);
I128_ASSERT_THAT(gcd(uint64_t{3}, int128_t{7}) == 1);
I128_ASSERT_THAT(gcd(uint64_t{3}, int128_t{-7}) == 1);
I128_ASSERT_THAT(gcd(uint64_t{3}, int128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(gcd(uint64_t{3}, int128_t{18446744073709551557ULL} * (-3)) == 3);
I128_ASSERT_THAT(gcd(uint64_t{3} * 2 * 5 * 7 * 11 * 13 * 17 * 19,
int128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(gcd(uint64_t{1000000000000000009ULL},
int128_t{1000000000000000009LL} * 1000000000000000009LL) ==
I128_ASSERT_THAT(math_functions::gcd(uint64_t{2}, int128_t{4}) == 2);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{2}, int128_t{-4}) == 2);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{3}, int128_t{7}) == 1);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{3}, int128_t{-7}) == 1);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{3}, int128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{3}, int128_t{18446744073709551557ULL} * (-3)) ==
3);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{3} * 2 * 5 * 7 * 11 * 13 * 17 * 19,
int128_t{18446744073709551557ULL} * 3) == 3);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{1000000000000000009ULL},
int128_t{1000000000000000009LL} * 1000000000000000009LL) ==
1000000000000000009ULL);
I128_ASSERT_THAT(gcd(uint64_t{0}, int128_t{1000000000000000009LL} * 1000000000000000009LL) ==
uint128_t{1000000000000000009LL} * 1000000000000000009ULL);
I128_ASSERT_THAT(gcd(uint64_t{18446744073709551557ULL}, int128_t{0}) ==
I128_ASSERT_THAT(
math_functions::gcd(uint64_t{0}, int128_t{1000000000000000009LL} * 1000000000000000009LL) ==
uint128_t{1000000000000000009LL} * 1000000000000000009ULL);
I128_ASSERT_THAT(math_functions::gcd(uint64_t{18446744073709551557ULL}, int128_t{0}) ==
18446744073709551557ULL);

// clang-format off

I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551557ULL}, int128_t{18446744073709551521ULL}) == int128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551533ULL}, int128_t{18446744073709551557ULL}) == int128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551521ULL}, int128_t{18446744073709551533ULL}) == int128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551557ULL}, int128_t{18446744073709551557ULL}) == int128_t{18446744073709551557ULL});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551533ULL}, int128_t{18446744073709551533ULL}) == int128_t{18446744073709551533ULL});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551521ULL}, int128_t{18446744073709551521ULL}) == int128_t{18446744073709551521ULL});

I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551557ULL}, int128_t{18446744073709551521ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551533ULL}, int128_t{18446744073709551557ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551521ULL}, int128_t{18446744073709551533ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551557ULL}, int128_t{18446744073709551557ULL}) == uint128_t{18446744073709551557ULL});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551533ULL}, int128_t{18446744073709551533ULL}) == uint128_t{18446744073709551533ULL});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551521ULL}, int128_t{18446744073709551521ULL}) == uint128_t{18446744073709551521ULL});

I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551557ULL}, uint128_t{18446744073709551521ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551533ULL}, uint128_t{18446744073709551557ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551521ULL}, uint128_t{18446744073709551533ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551557ULL}, uint128_t{18446744073709551557ULL}) == uint128_t{18446744073709551557ULL});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551533ULL}, uint128_t{18446744073709551533ULL}) == uint128_t{18446744073709551533ULL});
I128_ASSERT_THAT(math_functions::gcd(int128_t{18446744073709551521ULL}, uint128_t{18446744073709551521ULL}) == uint128_t{18446744073709551521ULL});

I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551557ULL}, uint128_t{18446744073709551521ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551533ULL}, uint128_t{18446744073709551557ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551521ULL}, uint128_t{18446744073709551533ULL}) == uint128_t{1});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551557ULL}, uint128_t{18446744073709551557ULL}) == uint128_t{18446744073709551557ULL});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551533ULL}, uint128_t{18446744073709551533ULL}) == uint128_t{18446744073709551533ULL});
I128_ASSERT_THAT(math_functions::gcd(uint128_t{18446744073709551521ULL}, uint128_t{18446744073709551521ULL}) == uint128_t{18446744073709551521ULL});

// clang-format on

ASSERT_THAT(math_functions::popcount(0U) == 0);
ASSERT_THAT(math_functions::popcount(1U << 1U) == 1);
ASSERT_THAT(math_functions::popcount(1U << 2U) == 1);
Expand Down

0 comments on commit 67a3cda

Please sign in to comment.