Skip to content

Commit

Permalink
Merge pull request #373 from ckormanyos/repair_gcd
Browse files Browse the repository at this point in the history 
Fix #374 via repair gcd
  • Loading branch information
@ckormanyos
ckormanyos authored Aug 14, 2023
2 parents 835daf2 + fe990b4 commit f32ee23
Show file tree
Hide file tree
Showing 4 changed files with 391 additions and 126 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ on how to use wide-integer.
- ![`example004_rootk_pow.cpp`](./examples/example004_rootk_pow.cpp) computes an integral seventh root and its corresponding power. A negative-valued cube root is also tested.
- ![`example005_powm.cpp`](./examples/example005_powm.cpp) tests the power-modulus function `powm`.
- ![`example005a_pow_factors_of_p99.cpp`](./examples/example005a_pow_factors_of_p99.cpp) verifies a beautiful, known prime factorization result from a classic tabulated value.
- ![`example006_gcd.cpp`](./examples/example006_gcd.cpp) tests the computation of a greatest common divisor `gcd`.
- ![`example006_gcd.cpp`](./examples/example006_gcd.cpp) tests several computations of greatest common divisor using the `gcd` function.
- ![`example007_random_generator.cpp`](./examples/example007_random_generator.cpp) computes some large pseudo-random integers.
- ![`example008_miller_rabin_prime.cpp`](./examples/example008_miller_rabin_prime.cpp) implements primality testing via Miller-Rabin.
- ![`example008a_miller_rabin_prime.cpp`](./examples/example008a_miller_rabin_prime.cpp) verifies Boost's interpretation of Miller-Rabin primality testing using `uintwide_t`-based types.
Expand Down
52 changes: 44 additions & 8 deletions examples/example006_gcd.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
///////////////////////////////////////////////////////////////////
// Copyright Christopher Kormanyos 2018 - 2022. //
// Copyright Christopher Kormanyos 2018 - 2023. //
// Distributed under the Boost Software License, //
// Version 1.0. (See accompanying file LICENSE_1_0.txt //
// or copy at http://www.boost.org/LICENSE_1_0.txt) //
Expand All @@ -20,16 +20,52 @@ auto ::math::wide_integer::example006_gcd() -> bool
using ::math::wide_integer::uint256_t;
#endif

WIDE_INTEGER_CONSTEXPR uint256_t a("0xD2690CD26CD57A3C443993851A70D3B62F841573668DF7B229508371A0AEDE7F");
WIDE_INTEGER_CONSTEXPR uint256_t b("0xFE719235CD0B1A314D4CA6940AEDC38BDF8E9484E68CE814EDAA17D87B0B4CC8");
auto result_is_ok = true;

WIDE_INTEGER_CONSTEXPR uint256_t c = gcd(a, b);
{
WIDE_INTEGER_CONSTEXPR auto a = uint256_t("0x1035452A5197CF882B5B5EB64C8CCFEE4D772F9F7B66A239649A43093464EFF5");
WIDE_INTEGER_CONSTEXPR auto b = uint256_t("0xB4F6151D727361113083D9A0DEB91B0B62A250F65DA6543823703D0140C873AD");

WIDE_INTEGER_CONSTEXPR auto result_is_ok = (static_cast<std::uint32_t>(c) == static_cast<std::uint32_t>(UINT32_C(12170749)));
WIDE_INTEGER_CONSTEXPR auto c = gcd(a, b);

#if defined(WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST) && (WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST != 0)
static_assert(result_is_ok, "Error: example006_gcd not OK!");
#endif
WIDE_INTEGER_CONSTEXPR auto result_gcd_is_ok = (static_cast<std::uint32_t>(c) == static_cast<std::uint32_t>(UINT32_C(11759761)));

#if defined(WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST) && (WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST != 0)
static_assert(result_gcd_is_ok, "Error: example006_gcd not OK!");
#endif

result_is_ok = (result_gcd_is_ok && result_is_ok);
}

{
WIDE_INTEGER_CONSTEXPR auto a = uint256_t("0xD2690CD26CD57A3C443993851A70D3B62F841573668DF7B229508371A0AEDE7F");
WIDE_INTEGER_CONSTEXPR auto b = uint256_t("0xFE719235CD0B1A314D4CA6940AEDC38BDF8E9484E68CE814EDAA17D87B0B4CC8");

WIDE_INTEGER_CONSTEXPR auto c = gcd(a, b);

WIDE_INTEGER_CONSTEXPR auto result_gcd_is_ok = (static_cast<std::uint32_t>(c) == static_cast<std::uint32_t>(UINT32_C(12170749)));

#if defined(WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST) && (WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST != 0)
static_assert(result_gcd_is_ok, "Error: example006_gcd not OK!");
#endif

result_is_ok = (result_gcd_is_ok && result_is_ok);
}

{
WIDE_INTEGER_CONSTEXPR auto a = uint256_t("0x7495AFF66DCB1085DC4CC294ECCBB1B455F65765DD4E9735564FDD80A05168A");
WIDE_INTEGER_CONSTEXPR auto b = uint256_t("0x7A0543EF0705942D09962172ED5038814AE6EDF8EED2FC6C52CF317D253BC81F");

WIDE_INTEGER_CONSTEXPR auto c = gcd(a, b);

WIDE_INTEGER_CONSTEXPR auto result_gcd_is_ok = (static_cast<std::uint32_t>(c) == static_cast<std::uint32_t>(UINT32_C(13520497)));

#if defined(WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST) && (WIDE_INTEGER_CONSTEXPR_IS_COMPILE_TIME_CONST != 0)
static_assert(result_gcd_is_ok, "Error: example006_gcd not OK!");
#endif

result_is_ok = (result_gcd_is_ok && result_is_ok);
}

return result_is_ok;
}
Expand Down
213 changes: 107 additions & 106 deletions math/wide_integer/uintwide_t.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,6 +181,99 @@
return dest;
}

template <class UnsignedIntegralType>
WIDE_INTEGER_CONSTEXPR auto clz_unsafe(UnsignedIntegralType v) noexcept -> std::enable_if_t<( (std::is_integral<UnsignedIntegralType>::value)
&& (std::is_unsigned<UnsignedIntegralType>::value)), unsigned>
{
using local_unsigned_integral_type = UnsignedIntegralType;

auto yy = static_cast<local_unsigned_integral_type>(UINT8_C(0));

auto nn = static_cast<unsigned>(std::numeric_limits<local_unsigned_integral_type>::digits);

auto cc = // NOLINT(altera-id-dependent-backward-branch)
static_cast<unsigned>
(
std::numeric_limits<local_unsigned_integral_type>::digits / static_cast<int>(INT8_C(2))
);

do
{
yy = static_cast<local_unsigned_integral_type>(v >> cc);

if(yy != static_cast<local_unsigned_integral_type>(UINT8_C(0)))
{
nn -= cc;

v = yy;
}

cc >>= static_cast<unsigned>(UINT8_C(1));
}
while(cc != static_cast<unsigned>(UINT8_C(0))); // NOLINT(altera-id-dependent-backward-branch)

return
static_cast<unsigned>
(
static_cast<unsigned>(nn) - static_cast<unsigned>(v)
);
}

template<typename UnsignedIntegralType>
WIDE_INTEGER_CONSTEXPR auto ctz_unsafe(const UnsignedIntegralType v) noexcept -> std::enable_if_t<( (std::is_integral<UnsignedIntegralType>::value)
&& (std::is_unsigned<UnsignedIntegralType>::value)), unsigned>
{
using local_unsigned_integral_type = UnsignedIntegralType;

constexpr auto local_digits = static_cast<unsigned>(std::numeric_limits<local_unsigned_integral_type>::digits);

const auto clz_mask =
static_cast<local_unsigned_integral_type>
(
static_cast<local_unsigned_integral_type>(~v)
& static_cast<local_unsigned_integral_type>(v - static_cast<local_unsigned_integral_type>(UINT8_C(1)))
);

return static_cast<unsigned>(local_digits - clz_unsafe(clz_mask));
}

template<typename UnsignedIntegralType>
WIDE_INTEGER_CONSTEXPR auto gcd_unsafe(UnsignedIntegralType u, UnsignedIntegralType v) -> std::enable_if_t<( (std::is_integral<UnsignedIntegralType>::value) // NOLINT(altera-id-dependent-backward-branch)
&& (std::is_unsigned<UnsignedIntegralType>::value)), UnsignedIntegralType>
{
using local_unsigned_integral_type = UnsignedIntegralType;

// Handle cases having (u != 0) and (v != 0).
if(u == static_cast<local_unsigned_integral_type>(UINT8_C(0))) { return v; }

if(v == static_cast<local_unsigned_integral_type>(UINT8_C(0))) { return u; }

// Shift the greatest power of 2 dividing both u and v.
const auto trz = static_cast<unsigned>(ctz_unsafe(u));

const auto shift_amount = (std::min)(trz, ctz_unsafe(v));

v >>= shift_amount;
u >>= trz;

do
{
// Reduce the GCD.

v >>= ctz_unsafe(v);

if(u > v)
{
std::swap(u, v);
}

v -= u;
}
while(v != static_cast<local_unsigned_integral_type>(UINT8_C(0))); // NOLINT(altera-id-dependent-backward-branch)

return static_cast<local_unsigned_integral_type>(u << shift_amount);
}

#if(__cplusplus >= 201703L)
} // namespace math::wide_integer::detail
#else
Expand Down Expand Up @@ -5884,68 +5977,9 @@
namespace detail {

template<typename UnsignedShortType>
WIDE_INTEGER_CONSTEXPR auto integer_gcd_reduce_short(UnsignedShortType u, UnsignedShortType v) -> UnsignedShortType
{
// This implementation of GCD reduction is based on an
// adaptation of existing code from Boost.Multiprecision.

for(;;)
{
if(u > v)
{
std::swap(u, v);
}

if(u == v)
{
break;
}

v -= u;
v >>= detail::lsb_helper(v);
}

return u;
}

template<typename UnsignedLargeType>
WIDE_INTEGER_CONSTEXPR auto integer_gcd_reduce_large(UnsignedLargeType u, UnsignedLargeType v) -> UnsignedLargeType
WIDE_INTEGER_CONSTEXPR auto integer_gcd_reduce(UnsignedShortType u, UnsignedShortType v) -> UnsignedShortType
{
// This implementation of GCD reduction is based on an
// adaptation of existing code from Boost.Multiprecision.

using local_ularge_type = UnsignedLargeType;
using local_ushort_type = typename detail::uint_type_helper<static_cast<size_t>(std::numeric_limits<local_ularge_type>::digits / 2)>::exact_unsigned_type;

for(;;)
{
if(u > v)
{
std::swap(u, v);
}

if(u == v)
{
break;
}

if(v <= static_cast<local_ularge_type>((std::numeric_limits<local_ushort_type>::max)()))
{
u = integer_gcd_reduce_short(static_cast<local_ushort_type>(v),
static_cast<local_ushort_type>(u));

break;
}

v -= u;

while(static_cast<std::uint_fast8_t>(static_cast<std::uint_fast8_t>(v) & static_cast<std::uint_fast8_t>(UINT8_C(1))) == static_cast<std::uint_fast8_t>(UINT8_C(0))) // NOLINT(hicpp-signed-bitwise,altera-id-dependent-backward-branch)
{
v >>= static_cast<unsigned>(UINT8_C(1));
}
}

return u;
return detail::gcd_unsafe(u, v);
}

} // namespace detail
Expand Down Expand Up @@ -5975,18 +6009,16 @@
using local_size_type = typename local_wide_integer_type::representation_type::size_type;

if(u == v)
{
{ // NOLINT(bugprone-branch-clone)
// This handles cases having (u = v) and also (u = v = 0).
result = u; // LCOV_EXCL_LINE
}

if((static_cast<local_ushort_type>(v) == static_cast<local_ushort_type>(UINT8_C(0))) && (v == 0U))
else if((static_cast<local_ushort_type>(v) == static_cast<local_ushort_type>(UINT8_C(0))) && (v == 0U))
{
// This handles cases having (v = 0) with (u != 0).
result = u; // LCOV_EXCL_LINE
}

if((static_cast<local_ushort_type>(u) == static_cast<local_ushort_type>(UINT8_C(0))) && (u == 0U))
else if((static_cast<local_ushort_type>(u) == static_cast<local_ushort_type>(UINT8_C(0))) && (u == 0U))
{
// This handles cases having (u = 0) with (v != 0).
result = v;
Expand Down Expand Up @@ -6028,7 +6060,7 @@
const auto vs = *v.crepresentation().cbegin();
const auto us = *u.crepresentation().cbegin();

u = detail::integer_gcd_reduce_short(vs, us);
u = detail::integer_gcd_reduce(vs, us);
}
else
{
Expand All @@ -6051,7 +6083,7 @@
const local_ularge_type v_large = detail::make_large(*v.crepresentation().cbegin(), my_v_hi);
const local_ularge_type u_large = detail::make_large(*u.crepresentation().cbegin(), my_u_hi);

u = detail::integer_gcd_reduce_large(v_large, u_large);
u = detail::integer_gcd_reduce(v_large, u_large);
}

break;
Expand All @@ -6071,38 +6103,7 @@
WIDE_INTEGER_CONSTEXPR auto gcd(const UnsignedShortType& u, const UnsignedShortType& v) -> std::enable_if_t<( (std::is_integral<UnsignedShortType>::value)
&& (std::is_unsigned<UnsignedShortType>::value)), UnsignedShortType>
{
using local_unsigned_short_type = UnsignedShortType;

auto result = local_unsigned_short_type { };

if(u > v)
{
result = gcd(v, u);
}

if(u == v)
{
// This handles cases having (u = v) and also (u = v = 0).
result = u;
}

if(v == static_cast<local_unsigned_short_type>(UINT8_C(0)))
{
// This handles cases having (v = 0) with (u != 0).
result = u;
}

if(u == static_cast<local_unsigned_short_type>(UINT8_C(0)))
{
// This handles cases having (u = 0) with (v != 0).
result = v;
}
else
{
result = detail::integer_gcd_reduce_short(u, v);
}

return result;
return detail::gcd_unsafe(u, v);
}

namespace detail {
Expand Down Expand Up @@ -6397,7 +6398,7 @@

generator_result_type value = generator_result_type();

auto it = result.representation().begin(); // NOLINT(llvm-qualified-auto,readability-qualified-auto)
auto it = result.representation().begin(); // NOLINT(llvm-qualified-auto,readability-qualified-auto,altera-id-dependent-backward-branch)

auto j = static_cast<unsigned_fast_type>(UINT8_C(0));

Expand Down Expand Up @@ -6562,7 +6563,7 @@

const auto m0 = static_cast<std::uint64_t>(np % pp0);

if((m0 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce_large(m0, pp0) != static_cast<std::uint64_t>(UINT8_C(1))))
if((m0 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce(m0, pp0) != static_cast<std::uint64_t>(UINT8_C(1))))
{
return false;
}
Expand All @@ -6576,7 +6577,7 @@

const auto m1 = static_cast<std::uint64_t>(np % pp1);

if((m1 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce_large(m1, pp1) != static_cast<std::uint64_t>(UINT8_C(1))))
if((m1 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce(m1, pp1) != static_cast<std::uint64_t>(UINT8_C(1))))
{
return false;
}
Expand All @@ -6590,7 +6591,7 @@

const auto m2 = static_cast<std::uint64_t>(np % pp2);

if((m2 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce_large(m2, pp2) != static_cast<std::uint64_t>(UINT8_C(1))))
if((m2 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce(m2, pp2) != static_cast<std::uint64_t>(UINT8_C(1))))
{
return false;
}
Expand All @@ -6604,7 +6605,7 @@

const auto m3 = static_cast<std::uint64_t>(np % pp3);

if((m3 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce_large(m3, pp3) != static_cast<std::uint64_t>(UINT8_C(1))))
if((m3 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce(m3, pp3) != static_cast<std::uint64_t>(UINT8_C(1))))
{
return false;
}
Expand All @@ -6618,7 +6619,7 @@

const auto m4 = static_cast<std::uint64_t>(np % pp4);

if((m4 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce_large(m4, pp4) != static_cast<std::uint64_t>(UINT8_C(1))))
if((m4 == static_cast<std::uint64_t>(UINT8_C(0))) || (detail::integer_gcd_reduce(m4, pp4) != static_cast<std::uint64_t>(UINT8_C(1))))
{
return false;
}
Expand Down
Loading

0 comments on commit f32ee23

Please sign in to comment.