Fix loss of precision in operator% and remove all individual (and duplicated) impls

include/boost/decimal/decimal128_t.hpp

@@ -41,6 +41,7 @@
 #include <boost/decimal/detail/components.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 
@@ -173,8 +174,6 @@ BOOST_DECIMAL_EXPORT class decimal128_t final
 
     friend constexpr auto d128_div_impl(const decimal128_t& lhs, const decimal128_t& rhs, decimal128_t& q, decimal128_t& r) noexcept -> void;
 
-    friend constexpr auto d128_mod_impl(const decimal128_t& lhs, const decimal128_t& rhs, const decimal128_t& q, decimal128_t& r) noexcept -> void;
-
     template <typename T>
     friend constexpr auto ilogb(T d) noexcept
         BOOST_DECIMAL_REQUIRES_RETURN(detail::is_decimal_floating_point_v, T, int);
@@ -1654,14 +1653,6 @@ constexpr auto d128_div_impl(const decimal128_t& lhs, const decimal128_t& rhs, d
 #  pragma warning(pop)
 #endif
 
-constexpr auto d128_mod_impl(const decimal128_t& lhs, const decimal128_t& rhs, const decimal128_t& q, decimal128_t& r) noexcept -> void
-{
-    constexpr decimal128_t zero {0, 0};
-
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator+(const decimal128_t& lhs, const decimal128_t& rhs) noexcept -> decimal128_t
 {
     #ifndef BOOST_DECIMAL_FAST_MATH
@@ -1968,7 +1959,7 @@ constexpr auto operator%(const decimal128_t& lhs, const decimal128_t& rhs) noexc
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        d128_mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/decimal32_t.hpp

@@ -42,6 +42,7 @@
 #include <boost/decimal/detail/chars_format.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 
@@ -158,7 +159,6 @@ BOOST_DECIMAL_EXPORT class decimal32_t final // NOLINT(cppcoreguidelines-special
     friend constexpr auto to_decimal(Decimal val) noexcept -> TargetType;
 
     friend constexpr auto div_impl(decimal32_t lhs, decimal32_t rhs, decimal32_t& q, decimal32_t& r) noexcept -> void;
-    friend constexpr auto mod_impl(decimal32_t lhs, decimal32_t rhs, const decimal32_t& q, decimal32_t& r) noexcept -> void;
 
     template <typename T>
     friend constexpr auto ilogb(T d) noexcept
@@ -2003,15 +2003,6 @@ constexpr auto div_impl(const decimal32_t lhs, const decimal32_t rhs, decimal32_
     q = detail::generic_div_impl<decimal32_t>(lhs_components, rhs_components);
 }
 
-constexpr auto mod_impl(const decimal32_t lhs, const decimal32_t rhs, const decimal32_t& q, decimal32_t& r) noexcept -> void
-{
-    constexpr decimal32_t zero {0, 0};
-
-    // https://en.cppreference.com/w/cpp/numeric/math/fmod
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator/(const decimal32_t lhs, const decimal32_t rhs) noexcept -> decimal32_t
 {
     decimal32_t q {};
@@ -2142,7 +2133,7 @@ constexpr auto operator%(const decimal32_t lhs, const decimal32_t rhs) noexcept
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/decimal64_t.hpp

@@ -44,6 +44,7 @@
 #include <boost/decimal/detail/to_chars_result.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 
@@ -184,8 +185,6 @@ BOOST_DECIMAL_EXPORT class decimal64_t final
 
     friend constexpr auto d64_div_impl(decimal64_t lhs, decimal64_t rhs, decimal64_t& q, decimal64_t& r) noexcept -> void;
 
-    friend constexpr auto d64_mod_impl(decimal64_t lhs, decimal64_t rhs, const decimal64_t& q, decimal64_t& r) noexcept -> void;
-
     template <typename T>
     friend constexpr auto ilogb(T d) noexcept
         BOOST_DECIMAL_REQUIRES_RETURN(detail::is_decimal_floating_point_v, T, int);
@@ -1633,15 +1632,6 @@ constexpr auto d64_div_impl(const decimal64_t lhs, const decimal64_t rhs, decima
     q = detail::d64_generic_div_impl<decimal64_t>(lhs_components, rhs.to_components(), sign);
 }
 
-constexpr auto d64_mod_impl(const decimal64_t lhs, const decimal64_t rhs, const decimal64_t& q, decimal64_t& r) noexcept -> void
-{
-    constexpr decimal64_t zero {0, 0};
-
-    // https://en.cppreference.com/w/cpp/numeric/math/fmod
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator+(const decimal64_t lhs, const decimal64_t rhs) noexcept -> decimal64_t
 {
     #ifndef BOOST_DECIMAL_FAST_MATH
@@ -1978,7 +1968,7 @@ constexpr auto operator%(const decimal64_t lhs, const decimal64_t rhs) noexcept
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        d64_mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/decimal_fast128_t.hpp

@@ -41,6 +41,7 @@
 #include <boost/decimal/detail/components.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 
@@ -1315,14 +1316,6 @@ constexpr auto d128f_div_impl(const decimal_fast128_t& lhs, const decimal_fast12
     q = decimal_fast128_t(static_cast<int128::uint128_t>(res_sig), res_exp, sign);
 }
 
-constexpr auto d128f_mod_impl(const decimal_fast128_t& lhs, const decimal_fast128_t& rhs, const decimal_fast128_t& q, decimal_fast128_t& r) -> void
-{
-    constexpr decimal_fast128_t zero {0, 0};
-
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator/(const decimal_fast128_t& lhs, const decimal_fast128_t& rhs) noexcept -> decimal_fast128_t
 {
     decimal_fast128_t q {};
@@ -1417,7 +1410,7 @@ constexpr auto operator%(const decimal_fast128_t& lhs, const decimal_fast128_t&
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        d128f_mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/decimal_fast32_t.hpp

@@ -42,6 +42,7 @@
 #include <boost/decimal/detail/chars_format.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 #include <limits>
@@ -127,8 +128,6 @@ BOOST_DECIMAL_EXPORT class alignas(4) decimal_fast32_t final
 
     friend constexpr auto div_impl(decimal_fast32_t lhs, decimal_fast32_t rhs, decimal_fast32_t& q, decimal_fast32_t& r) noexcept -> void;
 
-    friend constexpr auto mod_impl(decimal_fast32_t lhs, decimal_fast32_t rhs, const decimal_fast32_t& q, decimal_fast32_t& r) noexcept -> void;
-
     // Attempts conversion to integral type:
     // If this is nan sets errno to EINVAL and returns 0
     // If this is not representable sets errno to ERANGE and returns 0
@@ -1297,15 +1296,6 @@ constexpr auto div_impl(const decimal_fast32_t lhs, const decimal_fast32_t rhs,
     q = detail::generic_div_impl<decimal_fast32_t>(lhs, rhs);
 }
 
-constexpr auto mod_impl(const decimal_fast32_t lhs, const decimal_fast32_t rhs, const decimal_fast32_t& q, decimal_fast32_t& r) noexcept -> void
-{
-    constexpr decimal_fast32_t zero {0, 0};
-
-    // https://en.cppreference.com/w/cpp/numeric/math/fmod
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator/(const decimal_fast32_t lhs, const decimal_fast32_t rhs) noexcept -> decimal_fast32_t
 {
     decimal_fast32_t q {};
@@ -1399,7 +1389,7 @@ constexpr auto operator%(const decimal_fast32_t lhs, const decimal_fast32_t rhs)
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/decimal_fast64_t.hpp

@@ -42,6 +42,7 @@
 #include <boost/decimal/detail/chars_format.hpp>
 #include <boost/decimal/detail/construction_sign.hpp>
 #include <boost/decimal/detail/from_chars_impl.hpp>
+#include <boost/decimal/detail/mod_impl.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
 
@@ -1440,15 +1441,6 @@ constexpr auto d64_fast_div_impl(const decimal_fast64_t& lhs, const decimal_fast
     q = decimal_fast64_t{res_sig, res_exp, sign};
 }
 
-constexpr auto d64_fast_mod_impl(const decimal_fast64_t lhs, const decimal_fast64_t rhs, const decimal_fast64_t& q, decimal_fast64_t& r) noexcept -> void
-{
-    constexpr decimal_fast64_t zero {0, 0};
-
-    // https://en.cppreference.com/w/cpp/numeric/math/fmod
-    auto q_trunc {q > zero ? floor(q) : ceil(q)};
-    r = lhs - (q_trunc * rhs);
-}
-
 constexpr auto operator/(const decimal_fast64_t& lhs, const decimal_fast64_t& rhs) noexcept -> decimal_fast64_t
 {
     decimal_fast64_t q {};
@@ -1548,7 +1540,7 @@ constexpr auto operator%(const decimal_fast64_t lhs, const decimal_fast64_t rhs)
 
     if (BOOST_DECIMAL_LIKELY(!isnan(q)))
     {
-        d64_fast_mod_impl(lhs, rhs, q, r);
+        detail::generic_mod_impl(lhs, lhs.to_components(), rhs, rhs.to_components(), q, r);
     }
 
     return r;

include/boost/decimal/detail/mod_impl.hpp

@@ -0,0 +1,90 @@
+// Copyright 2023 - 2026 Matt Borland
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#ifndef BOOST_DECIMAL_DETAIL_MOD_IMPL_HPP
+#define BOOST_DECIMAL_DETAIL_MOD_IMPL_HPP
+
+#include <boost/decimal/detail/config.hpp>
+#include <boost/decimal/detail/power_tables.hpp>
+#include <boost/decimal/detail/int128.hpp>
+#include <boost/decimal/detail/u256.hpp>
+#include <boost/decimal/detail/promotion.hpp>
+#include <boost/decimal/detail/type_traits.hpp>
+#include <boost/decimal/detail/concepts.hpp>
+
+namespace boost {
+namespace decimal {
+namespace detail {
+
+template <typename T>
+constexpr auto get_pow_10(const std::size_t value) noexcept -> T
+{
+    static_assert(std::is_same<T, std::uint64_t>::value, "Should not be used");
+    return static_cast<T>(value);
+}
+
+template <>
+constexpr auto get_pow_10(const std::size_t value) noexcept -> std::uint64_t
+{
+    return pow10<std::uint64_t>(value);
+}
+
+template <>
+constexpr auto get_pow_10(const std::size_t value) noexcept -> int128::uint128_t
+{
+    return pow10_u128(value);
+}
+
+template <>
+constexpr auto get_pow_10(const std::size_t value) noexcept -> u256
+{
+    return pow10_256(value);
+}
+
+template <typename DecimalType, typename ComponentsType>
+constexpr auto generic_mod_impl(const DecimalType& lhs, const ComponentsType& lhs_components,
+                                const DecimalType& rhs, const ComponentsType& rhs_components,
+                                const DecimalType& q, DecimalType& r) noexcept
+    BOOST_DECIMAL_REQUIRES_TWO_RETURN(is_decimal_floating_point_v, DecimalType, is_decimal_floating_point_components_v, ComponentsType, void)
+{
+    using promoted_integer_type = std::conditional_t<decimal_val_v<DecimalType> < 64, std::uint64_t,
+                                     std::conditional_t<decimal_val_v<DecimalType> < 128, int128::uint128_t, u256>>;
+
+    const auto common_exp {std::min(lhs_components.exp, rhs_components.exp)};
+    const auto lhs_scaling {lhs_components.exp - common_exp};
+    const auto rhs_scaling {rhs_components.exp - common_exp};
+
+    // An approximation of the most digits we can hold without actually having to count the digits
+    constexpr auto max_scaling {std::numeric_limits<promoted_integer_type>::digits10 - std::numeric_limits<DecimalType>::digits10};
+
+    if (std::max(lhs_scaling, rhs_scaling) <= max_scaling)
+    {
+        BOOST_DECIMAL_ASSERT(lhs_scaling >= 0);
+        BOOST_DECIMAL_ASSERT(rhs_scaling >= 0);
+
+        promoted_integer_type scaled_lhs {lhs_components.sig};
+        promoted_integer_type scaled_rhs {rhs_components.sig};
+
+        scaled_lhs *= get_pow_10<promoted_integer_type>(static_cast<std::size_t>(lhs_scaling));
+        scaled_rhs *= get_pow_10<promoted_integer_type>(static_cast<std::size_t>(rhs_scaling));
+
+        const auto remainder_coeff {scaled_lhs % scaled_rhs};
+
+        r = DecimalType{remainder_coeff, common_exp, lhs_components.sign};
+    }
+    else
+    {
+        constexpr DecimalType zero {0, 0};
+
+        // https://en.cppreference.com/w/cpp/numeric/math/fmod
+        const auto q_trunc {q > zero ? floor(q) : ceil(q)};
+        r = lhs - (q_trunc * rhs);
+    }
+}
+
+} // namespace detail
+} // namespace decimal
+} // namespace boost
+
+#endif // BOOST_DECIMAL_DETAIL_MOD_IMPL_HPP

include/boost/decimal/detail/power_tables.hpp

@@ -225,6 +225,12 @@ constexpr auto pow10(const boost::int128::uint128_t n) noexcept -> boost::int128
     return impl::BOOST_DECIMAL_DETAIL_INT128_pow10[static_cast<std::size_t>(n.low)];
 }
 
+constexpr auto pow10_u128(const std::size_t n) noexcept -> boost::int128::uint128_t
+{
+    BOOST_DECIMAL_ASSERT(n < 39);
+    return impl::BOOST_DECIMAL_DETAIL_INT128_pow10[n];
+}
+
 #ifdef BOOST_DECIMAL_HAS_INT128
 
 template <>

include/boost/decimal/detail/type_traits.hpp

@@ -9,6 +9,7 @@
 #include <boost/decimal/fwd.hpp>
 #include <boost/decimal/detail/config.hpp>
 #include <boost/decimal/detail/u256.hpp>
+#include <boost/decimal/detail/components.hpp>
 #include "int128.hpp"
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
@@ -18,7 +19,6 @@
 namespace boost {
 namespace decimal {
 namespace detail {
-
 template <typename T>
 struct is_signed { static constexpr bool value = std::is_signed<T>::value; };
 
@@ -166,6 +166,24 @@ constexpr bool is_decimal_floating_point<T>::value;
 template <typename T>
 constexpr bool is_decimal_floating_point_v = is_decimal_floating_point<T>::value;
 
+template <typename T>
+struct is_decimal_floating_point_components { static constexpr bool value = false; };
+
+template <>
+struct is_decimal_floating_point_components<decimal32_t_components> { static constexpr bool value = true; };
+
+template <>
+struct is_decimal_floating_point_components<decimal64_t_components> { static constexpr bool value = true; };
+
+template <>
+struct is_decimal_floating_point_components<decimal128_t_components> { static constexpr bool value = true; };
+
+template <typename T>
+constexpr bool is_decimal_floating_point_components<T>::value;
+
+template <typename T>
+constexpr bool is_decimal_floating_point_components_v =  is_decimal_floating_point_components<T>::value;
+
 template <typename T>
 struct is_ieee_type { static constexpr bool value = false; };
 

include/boost/decimal/detail/u256.hpp

@@ -1280,7 +1280,7 @@ class numeric_limits<boost::decimal::detail::u256>
     static constexpr bool is_bounded = true;
     static constexpr bool is_modulo = true;
     static constexpr int digits = 256;
-    static constexpr int digits10 = 76;
+    static constexpr int digits10 = 78;
     static constexpr int max_digits10 = 0;
     static constexpr int radix = 2;
     static constexpr int min_exponent = 0;

test/random_decimal128_math.cpp

@@ -531,6 +531,12 @@ int main()
 
     test_spot_sub(-813150, -905406);
 
+    BOOST_TEST_EQ(decimal128_t{"Inf"} / decimal128_t{-1000}, -std::numeric_limits<decimal128_t>::infinity());
+
+    BOOST_TEST_EQ(decimal128_t{2} % decimal128_t{1}, decimal128_t{0});
+    BOOST_TEST_EQ(decimal128_t{"123456789"} / decimal128_t{"1.00000000000001"}, decimal128_t{"123456788.9999987654321100000123457"});
+    BOOST_TEST_EQ(decimal128_t{"123456789"} % decimal128_t{"1.00000000000001"}, decimal128_t{"0.99999876543212"});
+
     return boost::report_errors();
 }