diff --git a/include/boost/multiprecision/cpp_dec_float.hpp b/include/boost/multiprecision/cpp_dec_float.hpp index e04417158..e51ebd50c 100644 --- a/include/boost/multiprecision/cpp_dec_float.hpp +++ b/include/boost/multiprecision/cpp_dec_float.hpp @@ -458,9 +458,143 @@ class cpp_dec_float prec_elem = (std::min)(cpp_dec_float_elem_number, (std::max)(elems, static_cast(2))); } static cpp_dec_float pow2(boost::long_long_type i); + exponent_type order() const { const bool bo_order_is_zero = ((!(isfinite)()) || (data[0] == static_cast(0u))); + + return (bo_order_is_zero ? static_cast(0) : static_cast(exp + order_prefix())); + } + + template + void serialize(Archive& ar, const unsigned int /*version*/) + { + for (unsigned i = 0; i < data.size(); ++i) + ar& boost::make_nvp("digit", data[i]); + ar& boost::make_nvp("exponent", exp); + ar& boost::make_nvp("sign", neg); + ar& boost::make_nvp("class-type", fpclass); + ar& boost::make_nvp("precision", prec_elem); + } + + void eval_round_self() + { + const bool needs_rounding = (((isfinite)() == true) && (iszero() == false)); + + if(needs_rounding) + { + using local_size_type = typename array_type::size_type; + using local_limb_type = typename array_type::value_type; + + // Obtain the number of base-10 digits on the zero'th limb. + const std::int32_t digits_limb_0 = static_cast(order_prefix()); + + constexpr std::int32_t local_max_digits10 = + std::int32_t + ( + cpp_dec_float_digits10 + 4 + ); + + const std::int32_t digits_limb_1_to_n = local_max_digits10 - digits_limb_0; + + // Find the index of the element that contains the least-significant base-10 digit. + std::int32_t least_digit_idx = std::int32_t( (digits_limb_1_to_n / cpp_dec_float_elem_digits10) + + ((digits_limb_1_to_n % cpp_dec_float_elem_digits10) != 0)); + + // Set the index of the element that contains the rounding base-10 digit. + const std::int32_t round_digit_idx = + std::int32_t(((digits_limb_1_to_n % cpp_dec_float_elem_digits10) != 0) + ? least_digit_idx + : least_digit_idx + 1); + + // Find the base-10 order (position) of the least-significant base-10 digit. + const std::int32_t least_digit_pos = + (std::int32_t(digits_limb_1_to_n % cpp_dec_float_elem_digits10) != 0) + ? std::int32_t + ( + cpp_dec_float_elem_digits10 + - std::int32_t(digits_limb_1_to_n % cpp_dec_float_elem_digits10) + ) + : 0; + + // Find the base-10 order (position) of the rounding base-10 digit. + const std::int32_t round_digit_pos = + std::int32_t((least_digit_pos != 0) ? least_digit_pos - 1 : cpp_dec_float_elem_digits10 - 1); + + // Get the value of the rounding base-10 digit. + const std::uint8_t round_digit_value = + detail::digit_at_pos_in_limb + ( + data[local_size_type(round_digit_idx)], + unsigned(round_digit_pos) + ); + + const local_limb_type least_digit_p10 = detail::pow10_maker_as_runtime_value(std::uint32_t(least_digit_pos)); + + // Clear the lower base-10 digits of the rounded element. + data[local_size_type(least_digit_idx)] -= local_limb_type(data[local_size_type(least_digit_idx)] % least_digit_p10); + + // Clear the lower base-10 limbs. + if(local_size_type(least_digit_idx + 1) < data.size() - 1U) + { + std::fill(data.begin() + local_size_type(least_digit_idx + 1), data.end(), local_limb_type(0U)); + } + + // Perform round-to-nearest with no tie-breaking whatsoever. + if(round_digit_value >= 5U) + { + data[local_size_type(least_digit_idx)] += least_digit_p10; + + // There is a carry from rounding up. + std::uint_fast8_t carry_out = + ((data[local_size_type(least_digit_idx)] >= cpp_dec_float_elem_mask) + ? static_cast(1U) + : static_cast(0U)); + + // Propogate the carry into the limbs of higher significance as needed. + if(carry_out != 0U) + { + data[local_size_type(least_digit_idx)] -= cpp_dec_float_elem_mask; + + --least_digit_idx; + + for( ; least_digit_idx >= 0 && (carry_out != 0U); --least_digit_idx) + { + const local_limb_type tt = local_limb_type(data[local_size_type(least_digit_idx)] + local_limb_type(carry_out)); + + carry_out = ((tt >= cpp_dec_float_elem_mask) ? static_cast(1U) + : static_cast(0U)); + + data[local_size_type(least_digit_idx)] = + static_cast(tt - ((carry_out != 0U) ? cpp_dec_float_elem_mask + : static_cast(0U))); + } + + if((least_digit_idx < 0) && (carry_out != 0U)) + { + // In rare cases, propagation of the carry reaches the zero'th limb + // of highest significance, and we must shift the data, create a new limb + // with the carry value of 1 and adjust the exponent accordingly. + std::copy_backward(data.cbegin(), + data.cend() - 1, + data.end()); + + data[0U] = carry_out; + + exp = static_cast(exp + static_cast(cpp_dec_float_elem_digits10)); + } + } + } + } + } + + private: + static bool data_elem_is_non_zero_predicate(const std::uint32_t& d) { return (d != static_cast(0u)); } + static bool data_elem_is_non_nine_predicate(const std::uint32_t& d) { return (d != static_cast(cpp_dec_float::cpp_dec_float_elem_mask - 1)); } + static bool char_is_nonzero_predicate(const char& c) { return (c != static_cast('0')); } + + exponent_type order_prefix() const + { // // Binary search to find the order of the leading term: // @@ -506,25 +640,9 @@ class cpp_dec_float } } - return (bo_order_is_zero ? static_cast(0) : static_cast(exp + prefix)); - } - - template - void serialize(Archive& ar, const unsigned int /*version*/) - { - for (unsigned i = 0; i < data.size(); ++i) - ar& boost::make_nvp("digit", data[i]); - ar& boost::make_nvp("exponent", exp); - ar& boost::make_nvp("sign", neg); - ar& boost::make_nvp("class-type", fpclass); - ar& boost::make_nvp("precision", prec_elem); + return prefix; } - private: - static bool data_elem_is_non_zero_predicate(const std::uint32_t& d) { return (d != static_cast(0u)); } - static bool data_elem_is_non_nine_predicate(const std::uint32_t& d) { return (d != static_cast(cpp_dec_float::cpp_dec_float_elem_mask - 1)); } - static bool char_is_nonzero_predicate(const char& c) { return (c != static_cast('0')); } - void from_unsigned_long_long(const boost::ulong_long_type u); template ::extract_ else { // Extract the data into an boost::ulong_long_type value. - cpp_dec_float xn(extract_integer_part()); + cpp_dec_float xn(*this); + + xn.eval_round_self(); + if (xn.isneg()) xn.negate(); + xn = xn.extract_integer_part(); + val = static_cast(xn.data[0]); const std::int32_t imax = (std::min)(static_cast(static_cast(xn.exp) / cpp_dec_float_elem_digits10), static_cast(cpp_dec_float_elem_number - static_cast(1))); @@ -1686,35 +1809,40 @@ boost::ulong_long_type cpp_dec_float::extract { return static_cast(extract_signed_long_long()); } - - if (exp < static_cast(0)) + else if (exp < static_cast(0)) { return static_cast(0u); } - - const cpp_dec_float xn(extract_integer_part()); - - boost::ulong_long_type val; - - if (xn.compare(ulong_long_max()) > 0) - { - return (std::numeric_limits::max)(); - } else { - // Extract the data into an boost::ulong_long_type value. - val = static_cast(xn.data[0]); + cpp_dec_float xn(*this); - const std::int32_t imax = (std::min)(static_cast(static_cast(xn.exp) / cpp_dec_float_elem_digits10), static_cast(cpp_dec_float_elem_number - static_cast(1))); + xn.eval_round_self(); - for (std::int32_t i = static_cast(1); i <= imax; i++) + xn = xn.extract_integer_part(); + + boost::ulong_long_type val; + + if (xn.compare(ulong_long_max()) > 0) { - val *= static_cast(cpp_dec_float_elem_mask); - val += static_cast(xn.data[i]); + return (std::numeric_limits::max)(); + } + else + { + // Extract the data into an boost::ulong_long_type value. + val = static_cast(xn.data[0]); + + const std::int32_t imax = (std::min)(static_cast(static_cast(xn.exp) / cpp_dec_float_elem_digits10), static_cast(cpp_dec_float_elem_number - static_cast(1))); + + for (std::int32_t i = static_cast(1); i <= imax; i++) + { + val *= static_cast(cpp_dec_float_elem_mask); + val += static_cast(xn.data[i]); + } } - } - return val; + return val; + } } template @@ -3089,32 +3217,52 @@ template inline void eval_floor(cpp_dec_float& result, const cpp_dec_float& x) { result = x; - if (!(x.isfinite)() || x.isint()) + + if((x.isfinite)() == false) { if ((x.isnan)()) errno = EDOM; - return; } + else + { + result.eval_round_self(); - if (x.isneg()) - result -= cpp_dec_float::one(); - result = result.extract_integer_part(); + if(result.isint() == false) + { + if(result.isneg()) + { + result -= cpp_dec_float::one(); + } + + result = result.extract_integer_part(); + } + } } template inline void eval_ceil(cpp_dec_float& result, const cpp_dec_float& x) { result = x; - if (!(x.isfinite)() || x.isint()) + + if ((x.isfinite)() == false) { if ((x.isnan)()) errno = EDOM; - return; } + else + { + result.eval_round_self(); - if (!x.isneg()) - result += cpp_dec_float::one(); - result = result.extract_integer_part(); + if(result.isint() == false) + { + if(result.isneg() == false) + { + result += cpp_dec_float::one(); + } + + result = result.extract_integer_part(); + } + } } template diff --git a/include/boost/multiprecision/detail/tables.hpp b/include/boost/multiprecision/detail/tables.hpp index 2e2fa90f2..c8bdbfb11 100644 --- a/include/boost/multiprecision/detail/tables.hpp +++ b/include/boost/multiprecision/detail/tables.hpp @@ -69,12 +69,41 @@ struct a029750 } }; -constexpr std::uint32_t pow10_maker(std::uint32_t n) +inline constexpr std::uint32_t pow10_maker(std::uint32_t n) { // Make the constant power of 10^n. return ((n == UINT32_C(0)) ? UINT32_C(1) : pow10_maker(n - UINT32_C(1)) * UINT32_C(10)); } +inline std::uint32_t pow10_maker_as_runtime_value(const std::uint32_t n) +{ + using local_p10_table_type = std::array; + + constexpr local_p10_table_type local_p10_table = + {{ + UINT32_C(1), + UINT32_C(10), + UINT32_C(100), + UINT32_C(1000), + UINT32_C(10000), + UINT32_C(100000), + UINT32_C(1000000), + UINT32_C(10000000), + UINT32_C(100000000), + UINT32_C(1000000000) + }}; + + return ((n < std::uint32_t(std::tuple_size::value)) + ? local_p10_table[typename local_p10_table_type::size_type(n)] + : local_p10_table.back()); +} + +template +inline constexpr std::uint8_t digit_at_pos_in_limb(const LimbType u, const unsigned pos) +{ + return std::uint8_t(LimbType(u / static_cast(pow10_maker(pos))) % LimbType(10U)); +} + }}}} // namespace boost::multiprecision::backends::detail #endif // BOOST_MP_DETAIL_TABLES_HPP diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 8b95430e8..ac69b4e69 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1129,6 +1129,7 @@ test-suite misc : [ run git_issue_265.cpp : : : [ check-target-builds ../config//has_mpfr : gmp mpfr : no ] ] [ run git_issue_277.cpp ] [ run git_issue_313.cpp ] + [ run git_issue_368.cpp ] [ compile git_issue_98.cpp : [ check-target-builds ../config//has_float128 : TEST_FLOAT128 quadmath : ] [ check-target-builds ../config//has_gmp : TEST_GMP gmp : ] diff --git a/test/git_issue_368.cpp b/test/git_issue_368.cpp new file mode 100644 index 000000000..6c7461517 --- /dev/null +++ b/test/git_issue_368.cpp @@ -0,0 +1,128 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2021 Christopher Kormanyos. 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) + +#include +#include + +#include + +#include "test.hpp" + +namespace local +{ + std::minstd_rand0 eng(static_cast(std::clock())); + + std::uniform_int_distribution dist_sign + ( + 0U, + 1U + ); + + void test_original_issue_368() + { + using decimal = boost::multiprecision::number, + boost::multiprecision::et_off>; + + decimal a = decimal {"69000"} / decimal {"184"}; + decimal b = ceil(a); + + BOOST_CHECK_EQUAL(b, 375); + } + + template + void test_floor_rounding() + { + using decimal_type = DecimalType; + + bool result_is_ok = true; + + for(std::int32_t lo_index = 101U; lo_index < 1010U; lo_index += 7U) + { + for(std::int32_t hi_index = 10001U; hi_index < 100010U; hi_index += 17U) + { + const bool lo_is_neg = ((dist_sign(eng) % 2U) == 0U); + const bool hi_is_neg = ((dist_sign(eng) % 2U) == 0U); + + const std::int32_t lo = ((lo_is_neg == false) ? lo_index : -lo_index); + const std::int32_t hi = ((hi_is_neg == false) ? hi_index : -hi_index); + + const std::int32_t lo_hi = lo * hi; + + const decimal_type a = decimal_type { lo_hi } / decimal_type { lo }; + const decimal_type b = floor(a); + + result_is_ok &= (static_cast(a) == b); + + if(result_is_ok == false) + { + break; + } + } + } + + BOOST_CHECK_EQUAL(result_is_ok, true); + } + + template + void test_ceil_rounding() + { + using decimal_type = DecimalType; + + bool result_is_ok = true; + + for(std::int32_t lo_index = 101U; lo_index < 1010U; lo_index += 7U) + { + for(std::int32_t hi_index = 10001U; hi_index < 100010U; hi_index += 17U) + { + const bool lo_is_neg = ((dist_sign(eng) % 2U) == 0U); + const bool hi_is_neg = ((dist_sign(eng) % 2U) == 0U); + + const std::int32_t lo = ((lo_is_neg == false) ? lo_index : -lo_index); + const std::int32_t hi = ((hi_is_neg == false) ? hi_index : -hi_index); + + const std::int32_t lo_hi = lo * hi; + + const decimal_type a = decimal_type { lo_hi } / decimal_type { lo }; + const decimal_type b = ceil(a); + + result_is_ok &= (static_cast(a) == b); + + if(result_is_ok == false) + { + break; + } + } + } + + BOOST_CHECK_EQUAL(result_is_ok, true); + } +} + +int main(int, char**) +{ + // Test the original code from the issue + // (or a simplified form thereof) + local::test_original_issue_368(); + + // Test floor/ceil for 10 decimal digits. + { + using decimal = boost::multiprecision::number, + boost::multiprecision::et_off>; + + local::test_floor_rounding(); + local::test_ceil_rounding(); + } + + // Test floor/ceil for 12 decimal digits. + { + using decimal = boost::multiprecision::number, + boost::multiprecision::et_off>; + + local::test_floor_rounding(); + local::test_ceil_rounding(); + } + + return boost::report_errors(); +}