From 4d43e0f5816fc13f79937a4389b9f77c6acaef52 Mon Sep 17 00:00:00 2001 From: Shen-Ta Hsieh Date: Tue, 29 Sep 2026 14:25:01 +0800 Subject: [PATCH] Round a value below the smallest subnormal in each rounding mode - The constructors of decimal32_t, decimal64_t and decimal128_t had a special branch for a significand of one digit below the smallest exponent. It rounded a tie up, and it ignored the upward and the downward modes below the first step. The branch is removed. The general branch after it rounds these values correctly. - coefficient_rounding set the coefficient to zero when all its digits drop, and it did not round. It now rounds a zero with a sticky bit through fenv_round, thus a directed mode can give the smallest subnormal. A zero input still gives zero. - decimal128_t called coefficient_rounding only for a coefficient type of 128 bits or more. A narrower coefficient far below the smallest exponent read past the end of the pow10 table. This is undefined behavior, and it can stop the program with SIGFPE. It now takes the same path. - For a 128-bit coefficient with few digits, coefficient_rounding divided in the narrow significand type. A shift of more than 9 or 19 digits did not fit in that type, thus the power of ten lost its high bits. The quotient was wrong, or the divide by zero stopped the program. The narrow divide now also needs a shift that fits. - Add a test for the four cases. --- include/boost/decimal/decimal128_t.hpp | 59 ++----------- include/boost/decimal/decimal32_t.hpp | 53 +----------- include/boost/decimal/decimal64_t.hpp | 54 +----------- .../boost/decimal/detail/fenv_rounding.hpp | 12 ++- test/Jamfile | 1 + test/github_issue_1476.cpp | 85 +++++++++++++++++++ 6 files changed, 104 insertions(+), 160 deletions(-) create mode 100644 test/github_issue_1476.cpp diff --git a/include/boost/decimal/decimal128_t.hpp b/include/boost/decimal/decimal128_t.hpp index 7c12eeec1..15fffc805 100644 --- a/include/boost/decimal/decimal128_t.hpp +++ b/include/boost/decimal/decimal128_t.hpp @@ -769,6 +769,12 @@ BOOST_DECIMAL_CUDA_CONSTEXPR decimal128_t::decimal128_t(T1 coeff, T2 exp, const coeff_digits = detail::coefficient_rounding(coeff, exp, biased_exp, is_negative, detail::num_digits(coeff)); } } + else if (biased_exp < -(detail::precision_v - 1)) + { + // A narrow coefficient far below the range also needs this rounding, + // else the pow10 call below reads past the end of its table + coeff_digits = detail::coefficient_rounding(coeff, exp, biased_exp, is_negative, detail::num_digits(coeff)); + } constexpr int128::uint128_t zero {0, 0}; auto reduced_coeff {static_cast(coeff)}; @@ -813,58 +819,7 @@ BOOST_DECIMAL_CUDA_CONSTEXPR decimal128_t::decimal128_t(T1 coeff, T2 exp, const const auto exp_delta {biased_exp - static_cast(detail::d128_max_biased_exponent)}; const auto digit_delta {coeff_digits - exp_delta}; - if (biased_exp < 0 && coeff_digits == 1) - { - // This needs to be flushed to 0 or rounded to subnormal min - rounding_mode current_round_mode {_boost_decimal_global_rounding_mode}; - - #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION - - if (!BOOST_DECIMAL_IS_CONSTANT_EVALUATED(coeff)) - { - current_round_mode = _boost_decimal_global_runtime_rounding_mode; - } - - #endif - - bool round {false}; - if (biased_exp == -1) - { - switch (current_round_mode) - { - case rounding_mode::fe_dec_to_nearest_from_zero: - BOOST_DECIMAL_FALLTHROUGH - case rounding_mode::fe_dec_to_nearest: - if (reduced_coeff >= 5U) - { - round = true; - } - break; - case rounding_mode::fe_dec_upward: - if (!is_negative && reduced_coeff != 0U) - { - round = true; - } - break; - default: - round = false; - break; - } - } - - if (round) - { - // Subnormal min is just 1 - bits_ = UINT64_C(1); - } - else - { - bits_ = UINT64_C(0); - } - - bits_.high |= is_negative ? detail::d128_sign_mask : UINT64_C(0); - } - else if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision_v) + if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision_v) { // Same overflow-fold pattern as d32/d64: post-shift coeff is <= max_significand_v // and biased_exp lands in [0, max], so pack_in_range routes to direct_pack. diff --git a/include/boost/decimal/decimal32_t.hpp b/include/boost/decimal/decimal32_t.hpp index 7f1854381..39e94b67a 100644 --- a/include/boost/decimal/decimal32_t.hpp +++ b/include/boost/decimal/decimal32_t.hpp @@ -703,58 +703,7 @@ BOOST_DECIMAL_CUDA_CONSTEXPR decimal32_t::decimal32_t(T1 coeff, T2 exp, const de const auto exp_delta {biased_exp - static_cast(detail::d32_max_biased_exponent)}; const auto digit_delta {coeff_digits - exp_delta}; - if (biased_exp < 0 && coeff_digits == 1) - { - // This needs to be flushed to 0 or rounded to subnormal min - rounding_mode current_round_mode {_boost_decimal_global_rounding_mode}; - - #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION - - if (!BOOST_DECIMAL_IS_CONSTANT_EVALUATED(coeff)) - { - current_round_mode = _boost_decimal_global_runtime_rounding_mode; - } - - #endif - - bool round {false}; - if (biased_exp == -1) - { - switch (current_round_mode) - { - case rounding_mode::fe_dec_to_nearest_from_zero: - BOOST_DECIMAL_FALLTHROUGH - case rounding_mode::fe_dec_to_nearest: - if (reduced_coeff >= 5U) - { - round = true; - } - break; - case rounding_mode::fe_dec_upward: - if (!is_negative && reduced_coeff != 0) - { - round = true; - } - break; - default: - round = false; - break; - } - } - - if (round) - { - // Subnormal min is just 1 - bits_ = UINT32_C(1); - } - else - { - bits_ = UINT32_C(0); - } - - bits_ |= is_negative ? detail::d32_sign_mask : UINT32_C(0); - } - else if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision) + if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision) { // After the shift, coeff has coeff_digits+digit_delta <= precision digits // (so coeff <= max_significand_v) and biased_exp lands in [0, max] because diff --git a/include/boost/decimal/decimal64_t.hpp b/include/boost/decimal/decimal64_t.hpp index 5a029d4d9..c7a2c7ae4 100644 --- a/include/boost/decimal/decimal64_t.hpp +++ b/include/boost/decimal/decimal64_t.hpp @@ -743,59 +743,7 @@ BOOST_DECIMAL_CUDA_CONSTEXPR decimal64_t::decimal64_t(T1 coeff, T2 exp, const de const auto exp_delta {biased_exp - static_cast(detail::d64_max_biased_exponent)}; const auto digit_delta {coeff_digits - exp_delta}; - if (biased_exp < 0 && coeff_digits == 1) - { - // This needs to be flushed to 0 or rounded to subnormal min - // e.g. 7e-399 should not become 70e-398 but 7e-400 should become 0 - rounding_mode current_round_mode {_boost_decimal_global_rounding_mode}; - - #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION - - if (!BOOST_DECIMAL_IS_CONSTANT_EVALUATED(coeff)) - { - current_round_mode = _boost_decimal_global_runtime_rounding_mode; - } - - #endif - - bool round {false}; - if (biased_exp == -1) - { - switch (current_round_mode) - { - case rounding_mode::fe_dec_to_nearest_from_zero: - BOOST_DECIMAL_FALLTHROUGH - case rounding_mode::fe_dec_to_nearest: - if (reduced_coeff >= 5U) - { - round = true; - } - break; - case rounding_mode::fe_dec_upward: - if (!is_negative && reduced_coeff != 0) - { - round = true; - } - break; - default: - round = false; - break; - } - } - - if (round) - { - // Subnormal min is just 1 - bits_ = UINT64_C(1); - } - else - { - bits_ = UINT64_C(0); - } - - bits_ |= is_negative ? detail::d64_sign_mask : UINT64_C(0); - } - else if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision_v) + if (digit_delta > 0 && coeff_digits + digit_delta <= detail::precision_v) { // Coeff stays in range (<= max_significand_v) by the branch's digit budget, // and biased_exp lands in [0, max] by construction. pack_in_range hits diff --git a/include/boost/decimal/detail/fenv_rounding.hpp b/include/boost/decimal/detail/fenv_rounding.hpp index 0ba39102e..2e3dddece 100644 --- a/include/boost/decimal/detail/fenv_rounding.hpp +++ b/include/boost/decimal/detail/fenv_rounding.hpp @@ -717,8 +717,13 @@ BOOST_DECIMAL_CUDA_CONSTEXPR auto coefficient_rounding(T1& coeff, T2& exp, T3& b if (BOOST_DECIMAL_UNLIKELY(shift > std::numeric_limits::digits10)) { - // Bounds check for our tables in pow10 - coeff = 0; + // Bounds check for our tables in pow10. All digits drop, so round a zero with a sticky + // bit in the current mode. A directed mode then gives the smallest subnormal. + demoted_integer_type zero_coeff {0U}; + const auto removed_digits {detail::fenv_round(zero_coeff, sign, coeff != 0U)}; + coeff = static_cast(zero_coeff); + exp += removed_digits + shift; + biased_exp += removed_digits + shift; return 1; } @@ -745,7 +750,8 @@ BOOST_DECIMAL_CUDA_CONSTEXPR auto coefficient_rounding(T1& coeff, T2& exp, T3& b // so prefer the int compare over a 256-bit compare. // This is slightly more conservative for the narrow band of (digits10+1)-digit values that fit in the // demoted type by virtue of its high-bit slack which land in the wide-divmod branch. - if (coeff_digits <= std::numeric_limits::digits10) + if (coeff_digits <= std::numeric_limits::digits10 && + shift <= std::numeric_limits::digits10) { const auto smaller_coeff {static_cast(coeff)}; const auto smaller_pow10 {static_cast(shift_pow_ten)}; diff --git a/test/Jamfile b/test/Jamfile index 3bd2e8116..0cdc85e53 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -142,6 +142,7 @@ run github_issue_1459_toward_zero.cpp : : : off ; run github_issue_1467.cpp ; run github_issue_1473.cpp ; +run github_issue_1476.cpp ; run link_1.cpp link_2.cpp link_3.cpp ; run quick.cpp ; diff --git a/test/github_issue_1476.cpp b/test/github_issue_1476.cpp new file mode 100644 index 000000000..a0ee9a740 --- /dev/null +++ b/test/github_issue_1476.cpp @@ -0,0 +1,85 @@ +// Copyright 2026 Shen-Ta Hsieh +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include +#include +#include + +using namespace boost::decimal; + +// A value below the smallest subnormal rounds in the current mode: a one digit significand, a +// value whose digits all drop, and a coefficient type narrower or wider than the significand. +template +void test() +{ + using sig_type = typename T::significand_type; + using wide_type = decimal128_t::significand_type; + constexpr int etiny {std::numeric_limits::min_exponent10 - std::numeric_limits::digits10 + 1}; + const T zero {0}; + + // Half of the smallest step is a tie, and the even neighbour is zero + fesetround(rounding_mode::fe_dec_to_nearest); + BOOST_TEST_EQ(T(sig_type{5}, etiny - 1), zero); + BOOST_TEST_EQ(T(sig_type{5}, etiny - 1, true), -zero); + BOOST_TEST(signbit(T(sig_type{5}, etiny - 1, true))); + BOOST_TEST_EQ(T(sig_type{5}, etiny + 1) * T(sig_type{1}, -2), zero); + // An int coefficient far below: decimal128_t read past the pow10 table here + BOOST_TEST_EQ(T(11, etiny - 70), zero); + BOOST_TEST_EQ(T(10, etiny - 55), zero); + // A 128-bit coefficient far below: the divisor of 32 and 64-bit types lost its high bits + BOOST_TEST_EQ(T(wide_type{123}, etiny - 33), zero); + BOOST_TEST_EQ(T(wide_type{UINT64_C(3976006815679288586)}, etiny - 38), zero); + + #ifndef BOOST_DECIMAL_NO_CONSTEVAL_DETECTION + // fesetround has an effect only with the detection of constant evaluation + const T smallest {std::numeric_limits::denorm_min()}; + + // A positive value goes up to the smallest step, also when all its digits drop + fesetround(rounding_mode::fe_dec_upward); + BOOST_TEST_EQ(T(sig_type{1}, etiny - 2), smallest); + BOOST_TEST_EQ(T(sig_type{1}, etiny + 1) * T(sig_type{1}, -3), smallest); + BOOST_TEST_EQ(T(sig_type{11}, etiny - 70), smallest); + BOOST_TEST_EQ(T(11, etiny - 70), smallest); + BOOST_TEST_EQ(T(wide_type{123}, etiny - 33), smallest); + // A zero stays zero + BOOST_TEST_EQ(T(sig_type{0}, etiny - 70), zero); + BOOST_TEST_EQ(T(sig_type{5}, etiny + 5) * T(sig_type{0}, -70), zero); + + // A negative value goes down to minus the smallest step, also when all its digits drop + fesetround(rounding_mode::fe_dec_downward); + BOOST_TEST_EQ(T(sig_type{1}, etiny - 2, true), -smallest); + BOOST_TEST_EQ(T(sig_type{7}, etiny - 1, true), -smallest); + BOOST_TEST_EQ(T(sig_type{1}, etiny + 1, true) * T(sig_type{1}, -3), -smallest); + BOOST_TEST_EQ(T(sig_type{11}, etiny - 70, true), -smallest); + + // A tie goes away from zero, but a value less than half of the smallest step goes to zero + fesetround(rounding_mode::fe_dec_to_nearest_from_zero); + BOOST_TEST_EQ(T(sig_type{5}, etiny - 1), smallest); + BOOST_TEST_EQ(T(sig_type{5}, etiny - 1, true), -smallest); + BOOST_TEST_EQ(T(sig_type{4}, etiny - 1), zero); + BOOST_TEST_EQ(T(sig_type{5}, etiny + 1) * T(sig_type{1}, -2), smallest); + BOOST_TEST_EQ(T(11, etiny - 70), zero); + BOOST_TEST_EQ(T(wide_type{123}, etiny - 33), zero); + + // All values go to zero and keep their sign + fesetround(rounding_mode::fe_dec_toward_zero); + BOOST_TEST_EQ(T(sig_type{9}, etiny - 1), zero); + BOOST_TEST(signbit(T(sig_type{9}, etiny - 1, true))); + BOOST_TEST_EQ(T(sig_type{1}, etiny + 1) * T(sig_type{9}, -2), zero); + BOOST_TEST_EQ(T(sig_type{11}, etiny - 70, true), -zero); + BOOST_TEST_EQ(T(-11, etiny - 70), -zero); + BOOST_TEST_EQ(T(wide_type{123}, etiny - 33), zero); + + fesetround(rounding_mode::fe_dec_to_nearest); + #endif +} + +int main() +{ + test(); + test(); + test(); + + return boost::report_errors(); +}