From e06fc7e19a164f318324f7b3266ecec85418a3f1 Mon Sep 17 00:00:00 2001 From: Juan Cruz Viotti Date: Fri, 7 Aug 2026 18:24:01 -0300 Subject: [PATCH 1/4] Fix Decimal division hang and make additive arithmetic exact Signed-off-by: Juan Cruz Viotti --- src/lang/numeric/big_coefficient.h | 140 +++++++++++++++----- src/lang/numeric/decimal.cc | 84 ++++++------ test/json/json_value_test.cc | 11 +- test/numeric/numeric_decimal_test.cc | 191 ++++++++++++++++++++++++++- 4 files changed, 338 insertions(+), 88 deletions(-) diff --git a/src/lang/numeric/big_coefficient.h b/src/lang/numeric/big_coefficient.h index 1597427d2..df70905b7 100644 --- a/src/lang/numeric/big_coefficient.h +++ b/src/lang/numeric/big_coefficient.h @@ -393,47 +393,119 @@ class BigCoefficient { return {std::move(quotient), std::move(remainder_big)}; } - auto remainder = this->clone(); - BigCoefficient quotient{this->length}; - quotient.length = this->length; - std::fill(quotient.words, quotient.words + quotient.length, 0ULL); - - while (remainder.compare(divisor) >= 0) { - auto remainder_top = static_cast( - remainder.words[remainder.length - 1]); - if (remainder.length > divisor.length) { - auto shift = remainder.length - divisor.length; - auto divisor_top = divisor.words[divisor.length - 1]; - auto estimate = - static_cast(remainder_top / (divisor_top + 1)); - if (estimate == 0) { - estimate = 1; + // Long division per Knuth, The Art of Computer Programming, volume 2, + // section 4.3.1, algorithm D, which normalizes both operands so that the + // top divisor word is at least half the base. This guarantees that each + // trial quotient word is at most one in excess after the two word + // correction test, so the cost is quadratic in the word count instead of + // linear in the magnitude of the quotient + auto normalizer = BASE / (divisor.words[divisor.length - 1] + 1); + + BigCoefficient normalized_divisor{divisor.length}; + normalized_divisor.length = divisor.length; + sourcemeta::core::uint128_t normalize_carry = 0; + for (std::uint32_t index = 0; index < divisor.length; index++) { + auto product = + static_cast(divisor.words[index]) * + normalizer + + normalize_carry; + normalized_divisor.words[index] = + static_cast(product % BASE); + normalize_carry = product / BASE; + } + + BigCoefficient normalized_dividend{this->length + 1}; + normalized_dividend.length = this->length + 1; + normalize_carry = 0; + for (std::uint32_t index = 0; index < this->length; index++) { + auto product = + static_cast(this->words[index]) * + normalizer + + normalize_carry; + normalized_dividend.words[index] = + static_cast(product % BASE); + normalize_carry = product / BASE; + } + + normalized_dividend.words[this->length] = + static_cast(normalize_carry); + + auto quotient_length = this->length - divisor.length + 1; + BigCoefficient quotient{quotient_length}; + quotient.length = quotient_length; + + auto top_divisor_word = normalized_divisor.words[divisor.length - 1]; + auto next_divisor_word = normalized_divisor.words[divisor.length - 2]; + + for (auto position = quotient_length; position > 0;) { + position--; + auto numerator = + static_cast( + normalized_dividend.words[position + divisor.length]) * + BASE + + normalized_dividend.words[position + divisor.length - 1]; + auto trial_word = numerator / top_divisor_word; + auto trial_remainder = numerator % top_divisor_word; + while (trial_word >= BASE || + trial_word * next_divisor_word > + trial_remainder * BASE + + normalized_dividend.words[position + divisor.length - 2]) { + trial_word--; + trial_remainder += top_divisor_word; + if (trial_remainder >= BASE) { + break; } + } - BigCoefficient estimate_big{1}; - estimate_big.words[0] = estimate; - estimate_big.length = 1; - - auto scaled = estimate_big.multiply_pow10(shift * BASE_DIGITS); - auto product = scaled.multiply(divisor); - - if (product.compare(remainder) > 0) { - remainder = remainder.subtract(divisor); - quotient.words[0]++; + std::uint64_t borrow = 0; + for (std::uint32_t index = 0; index < divisor.length; index++) { + auto subtrahend = trial_word * normalized_divisor.words[index] + borrow; + auto subtrahend_low = static_cast(subtrahend % BASE); + borrow = static_cast(subtrahend / BASE); + auto &word = normalized_dividend.words[position + index]; + if (word < subtrahend_low) { + word += BASE - subtrahend_low; + borrow++; } else { - remainder = remainder.subtract(product); - BigCoefficient estimated_quotient{shift + 1}; - std::fill(estimated_quotient.words, estimated_quotient.words + shift, - 0ULL); - estimated_quotient.words[shift] = estimate; - estimated_quotient.length = shift + 1; - quotient = quotient.add(estimated_quotient); + word -= subtrahend_low; } + } + auto top_word = normalized_dividend.words[position + divisor.length]; + if (top_word < borrow) { + trial_word--; + std::uint64_t add_carry = 0; + for (std::uint32_t index = 0; index < divisor.length; index++) { + auto &word = normalized_dividend.words[position + index]; + auto sum = word + normalized_divisor.words[index] + add_carry; + if (sum >= BASE) { + word = sum - BASE; + add_carry = 1; + } else { + word = sum; + add_carry = 0; + } + } + + normalized_dividend.words[position + divisor.length] = + top_word + add_carry - borrow; } else { - remainder = remainder.subtract(divisor); - quotient.words[0]++; + normalized_dividend.words[position + divisor.length] = + top_word - borrow; } + + quotient.words[position] = static_cast(trial_word); + } + + BigCoefficient remainder{divisor.length}; + remainder.length = divisor.length; + sourcemeta::core::uint128_t denormalize_carry = 0; + for (auto index = divisor.length; index > 0; index--) { + auto current = + denormalize_carry * BASE + normalized_dividend.words[index - 1]; + remainder.words[index - 1] = + static_cast(current / normalizer); + denormalize_carry = current % normalizer; } quotient.trim(); diff --git a/src/lang/numeric/decimal.cc b/src/lang/numeric/decimal.cc index d0ebb5d94..765f1192a 100644 --- a/src/lang/numeric/decimal.cc +++ b/src/lang/numeric/decimal.cc @@ -877,6 +877,13 @@ auto Decimal::to_uint32() const -> std::uint32_t { } auto Decimal::to_float() const -> float { + // IEEE 754-2019 section 6.2 requires an operation that signals an invalid + // operation exception and delivers a floating point result to deliver a + // quiet NaN, which is the case when converting a signaling NaN + if (this->is_nan()) { + return std::numeric_limits::quiet_NaN(); + } + try { return std::stof(this->to_scientific_string()); } catch (const std::out_of_range &) { @@ -885,6 +892,13 @@ auto Decimal::to_float() const -> float { } auto Decimal::to_double() const -> double { + // IEEE 754-2019 section 6.2 requires an operation that signals an invalid + // operation exception and delivers a floating point result to deliver a + // quiet NaN, which is the case when converting a signaling NaN + if (this->is_nan()) { + return std::numeric_limits::quiet_NaN(); + } + try { return std::stod(this->to_scientific_string()); } catch (const std::out_of_range &) { @@ -1791,8 +1805,6 @@ auto Decimal::operator+=(const Decimal &other) -> Decimal & { store_big_result(this->coefficient_, this->coefficient_high_, this->flags_, std::move(result_big), result_negative); this->exponent_ = result_exponent; - round_to_precision(this->coefficient_, this->coefficient_high_, - this->exponent_, this->flags_); return *this; } @@ -1993,52 +2005,38 @@ auto Decimal::operator%=(const Decimal &other) -> Decimal & { return *this; } - Decimal quotient{*this}; - quotient /= other; + // The General Decimal Arithmetic Specification defines remainder as "the + // residue of the dividend after the operation of calculating integer + // division" and states that "the sign of the result, if non-zero, is the + // same as that of the original dividend", so the result is derived from + // the exact big integer division rather than from rounded arithmetic + Decimal dividend_magnitude{*this}; + dividend_magnitude.flags_ = + static_cast(dividend_magnitude.flags_ & ~FLAG_SIGN); + Decimal divisor_magnitude{other}; + divisor_magnitude.flags_ = + static_cast(divisor_magnitude.flags_ & ~FLAG_SIGN); + if (dividend_magnitude < divisor_magnitude) { + return *this; + } - if (quotient.is_finite() && !quotient.is_zero()) { - if (quotient.exponent_ < 0) { - if (quotient.flags_ & FLAG_BIG) { - auto digit_string = coefficient_to_digit_string( - quotient.coefficient_, quotient.coefficient_high_, quotient.flags_); - auto number_of_digits = static_cast(digit_string.size()); - auto digits_to_remove = -quotient.exponent_; - if (digits_to_remove >= number_of_digits) { - quotient = Decimal{}; - } else { - auto integer_string = digit_string.substr( - 0, static_cast(number_of_digits - digits_to_remove)); - auto old_sign = - static_cast(quotient.flags_ & FLAG_SIGN); - // The assignment below releases the current coefficient, so freeing - // it explicitly here as well would free the same allocation twice - quotient = Decimal{integer_string}; - quotient.flags_ = - static_cast(quotient.flags_ | old_sign); - } + bool result_negative = (this->flags_ & FLAG_SIGN) != 0; + auto result_exponent = std::min(this->exponent_, other.exponent_); - } else { - auto coefficient = quotient.coefficient_; - auto exponent = quotient.exponent_; - while (exponent < 0 && coefficient > 0) { - coefficient /= 10; - exponent++; - } + auto dividend_big = coefficient_as_big(this->coefficient_, + this->coefficient_high_, this->flags_); + auto divisor_big = coefficient_as_big(other.coefficient_, + other.coefficient_high_, other.flags_); - if (exponent < 0) { - quotient = Decimal{}; - } else { - quotient.coefficient_ = coefficient; - quotient.exponent_ = exponent; - } - } - } - } + BigCoefficient::align_exponents(dividend_big, divisor_big, this->exponent_, + other.exponent_); - Decimal product{quotient}; - product *= other; - *this -= product; + auto [quotient, remainder] = dividend_big.divide_modulo(divisor_big); + free_big_coefficient(this->coefficient_, this->flags_); + store_big_result(this->coefficient_, this->coefficient_high_, this->flags_, + std::move(remainder), result_negative); + this->exponent_ = result_exponent; return *this; } diff --git a/test/json/json_value_test.cc b/test/json/json_value_test.cc index bcde80a37..1179e8dc8 100644 --- a/test/json/json_value_test.cc +++ b/test/json/json_value_test.cc @@ -847,24 +847,17 @@ TEST(deep_copy_of_a_nested_object) { TEST(add_integer_overflow_promotes_to_decimal) { const sourcemeta::core::JSON left{std::numeric_limits::min()}; const sourcemeta::core::JSON right{-1}; - // The integer sum overflows int64, so the operator promotes to Decimal. The - // exact value would be -9223372036854775809, but the Decimal add currently - // rounds to the working precision (tracked in the decimal precision bug - // report). Pin the exact current output so a change in either direction is - // caught const auto result{left + right}; EXPECT_TRUE(result.is_decimal()); - EXPECT_EQ(result.to_decimal().to_string(), "-9.223372036854776e+18"); + EXPECT_EQ(result.to_decimal().to_string(), "-9223372036854775809"); } TEST(subtract_integer_overflow_promotes_to_decimal) { const sourcemeta::core::JSON left{std::numeric_limits::max()}; const sourcemeta::core::JSON right{-1}; - // See the note on the addition overflow test above. Exact value would be - // 9223372036854775808; pin the current rounded output const auto result{left - right}; EXPECT_TRUE(result.is_decimal()); - EXPECT_EQ(result.to_decimal().to_string(), "9.223372036854776e+18"); + EXPECT_EQ(result.to_decimal().to_string(), "9223372036854775808"); } TEST(copy_self_assignment) { diff --git a/test/numeric/numeric_decimal_test.cc b/test/numeric/numeric_decimal_test.cc index b3e6ea3d2..4900cb866 100644 --- a/test/numeric/numeric_decimal_test.cc +++ b/test/numeric/numeric_decimal_test.cc @@ -343,7 +343,7 @@ TEST(large_multiplication) { TEST(high_precision_addition) { const sourcemeta::core::Decimal left{"0.1111111111111111111111111111"}; const sourcemeta::core::Decimal right{"0.2222222222222222222222222222"}; - const sourcemeta::core::Decimal expected{"0.3333333333333333"}; + const sourcemeta::core::Decimal expected{"0.3333333333333333333333333333"}; EXPECT_EQ(left + right, expected); } @@ -813,6 +813,18 @@ TEST(to_double_nan) { EXPECT_TRUE(std::isnan(result)); } +TEST(to_float_signaling_nan) { + const auto nan_value{sourcemeta::core::Decimal::snan()}; + const float result{nan_value.to_float()}; + EXPECT_TRUE(std::isnan(result)); +} + +TEST(to_double_signaling_nan) { + const auto nan_value{sourcemeta::core::Decimal::snan()}; + const double result{nan_value.to_double()}; + EXPECT_TRUE(std::isnan(result)); +} + TEST(to_double_infinity) { const auto inf_value{sourcemeta::core::Decimal::infinity()}; const double result{inf_value.to_double()}; @@ -4126,7 +4138,7 @@ TEST(divide_integer_negative_big_quotient) { TEST(subtract_with_borrow_cascade) { const sourcemeta::core::Decimal left{"1000000000000000000000000000000"}; const sourcemeta::core::Decimal right{1}; - EXPECT_EQ((left - right).to_string(), "1.0000000000000000e+30"); + EXPECT_EQ((left - right).to_string(), "999999999999999999999999999999"); } TEST(subtract_bigs_to_small_result) { @@ -4295,3 +4307,178 @@ TEST(divide_integer_shifted_estimate_small_quotient) { "1000000000000000000000000000000000000000000"}; EXPECT_EQ(dividend.divide_integer(divisor).to_string(), "999"); } + +TEST(divide_integer_multi_word_divisor_large_quotient) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ(dividend.divide_integer(divisor).to_string(), + "810000007290000066339"); +} + +TEST(divide_integer_multi_word_divisor_exact_multiple) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999998704799994976289954955"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ(dividend.divide_integer(divisor).to_string(), + "810000007290000066339"); +} + +TEST(divide_integer_two_word_divisor_power_of_ten_quotient) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal divisor{ + "999999999999999999000000000000000005"}; + EXPECT_EQ(dividend.divide_integer(divisor).to_string(), "1000000000"); +} + +TEST(divide_integer_heap_dividend_multi_word_divisor) { + const sourcemeta::core::Decimal dividend{ + "1428571428571428571428571428571428571428571428571428571428571428571428" + "57142857142857142857142857142"}; + const sourcemeta::core::Decimal divisor{ + "9999999999999999999999999999999999999997"}; + EXPECT_EQ(dividend.divide_integer(divisor).to_string(), + "14285714285714285714285714285714285714290000000000000000000"); +} + +TEST(divide_multi_word_divisor_large_quotient) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ(dividend / divisor, + sourcemeta::core::Decimal{"8.100000072900001e+20"}); +} + +TEST(divisible_by_multi_word_divisor_false) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_FALSE(dividend.divisible_by(divisor)); +} + +TEST(divisible_by_multi_word_divisor_true) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999998704799994976289954955"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_TRUE(dividend.divisible_by(divisor)); +} + +// The trial quotient word passes the two word correction test because the +// middle divisor word is zero, yet is still one in excess due to the lowest +// divisor word, forcing the rare add back step of the long division +TEST(divide_integer_trial_word_excess_add_back) { + const sourcemeta::core::Decimal dividend{ + "499999999999999999500000000000000000000000000000000000" + "999999999999999998"}; + const sourcemeta::core::Decimal divisor{ + "500000000000000000000000000000000000000000000000000001"}; + EXPECT_EQ(dividend.divide_integer(divisor).to_string(), "999999999999999998"); +} + +TEST(modulo_trial_word_excess_add_back) { + const sourcemeta::core::Decimal dividend{ + "499999999999999999500000000000000000000000000000000000" + "999999999999999998"}; + const sourcemeta::core::Decimal divisor{ + "500000000000000000000000000000000000000000000000000001"}; + EXPECT_EQ((dividend % divisor).to_string(), + "500000000000000000000000000000000000000000000000000000"); +} + +TEST(modulo_eighteen_digit_exact_multiple) { + const sourcemeta::core::Decimal dividend{"999999999999999999"}; + const sourcemeta::core::Decimal divisor{7}; + EXPECT_EQ((dividend % divisor).to_string(), "0"); +} + +TEST(modulo_seventeen_digit_dividend) { + const sourcemeta::core::Decimal dividend{"99999999999999999"}; + const sourcemeta::core::Decimal divisor{7}; + EXPECT_EQ((dividend % divisor).to_string(), "4"); +} + +TEST(modulo_negative_seventeen_digit_dividend) { + const sourcemeta::core::Decimal dividend{"-99999999999999999"}; + const sourcemeta::core::Decimal divisor{7}; + EXPECT_EQ((dividend % divisor).to_string(), "-4"); +} + +TEST(modulo_exact_multiple_takes_dividend_sign) { + const sourcemeta::core::Decimal dividend{"-999999999999999999"}; + const sourcemeta::core::Decimal divisor{7}; + const sourcemeta::core::Decimal result{dividend % divisor}; + EXPECT_TRUE(result.is_zero()); + EXPECT_TRUE(result.is_signed()); +} + +TEST(modulo_multi_word_divisor) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ((dividend % divisor).to_string(), "1295200005023710045044"); +} + +TEST(modulo_multi_word_divisor_exact_multiple) { + const sourcemeta::core::Decimal dividend{ + "999999999999999999999998704799994976289954955"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ((dividend % divisor).to_string(), "0"); +} + +TEST(modulo_heap_dividend_multi_word_divisor) { + const sourcemeta::core::Decimal dividend{ + "1428571428571428571428571428571428571428571428571428571428571428571428" + "57142857142857142857142857142"}; + const sourcemeta::core::Decimal divisor{ + "9999999999999999999999999999999999999997"}; + EXPECT_EQ((dividend % divisor).to_string(), "12857142857142857142"); +} + +TEST(modulo_high_precision_fractional_operands) { + const sourcemeta::core::Decimal dividend{"1.000000000000000001"}; + const sourcemeta::core::Decimal divisor{"0.000000000000000002"}; + EXPECT_EQ(dividend % divisor, + sourcemeta::core::Decimal{"0.000000000000000001"}); +} + +TEST(modulo_huge_divisor_returns_dividend) { + const sourcemeta::core::Decimal dividend{"0.5"}; + const sourcemeta::core::Decimal divisor{"1e+100"}; + EXPECT_EQ(dividend % divisor, sourcemeta::core::Decimal{"0.5"}); +} + +TEST(add_exact_below_int64_minimum) { + const sourcemeta::core::Decimal left{ + std::numeric_limits::min()}; + const sourcemeta::core::Decimal right{-1}; + EXPECT_EQ((left + right).to_string(), "-9223372036854775809"); +} + +TEST(subtract_exact_above_int64_maximum) { + const sourcemeta::core::Decimal left{ + std::numeric_limits::max()}; + const sourcemeta::core::Decimal right{-1}; + EXPECT_EQ((left - right).to_string(), "9223372036854775808"); +} + +TEST(subtract_exact_twenty_five_digit_result) { + const sourcemeta::core::Decimal left{"10000000000000000000000000"}; + const sourcemeta::core::Decimal right{1}; + EXPECT_EQ((left - right).to_string(), "9999999999999999999999999"); +} + +TEST(add_exact_heap_operands) { + const sourcemeta::core::Decimal left{ + "999999999999999999999999999999999999999999999"}; + const sourcemeta::core::Decimal right{ + "123456789012345678901234567890123456789012345"}; + EXPECT_EQ((left + right).to_string(), + "1123456789012345678901234567890123456789012344"); +} + +TEST(add_exact_fractional_carry_chain) { + const sourcemeta::core::Decimal left{"9999999999999999.9999999999999999"}; + const sourcemeta::core::Decimal right{"0.0000000000000001"}; + EXPECT_EQ(left + right, sourcemeta::core::Decimal{"10000000000000000"}); +} From d6efc05930d898665e2880d38f22cc0260b64eac Mon Sep 17 00:00:00 2001 From: Juan Cruz Viotti Date: Fri, 7 Aug 2026 19:05:52 -0300 Subject: [PATCH 2/4] Fix Signed-off-by: Juan Cruz Viotti --- src/lang/numeric/big_coefficient.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lang/numeric/big_coefficient.h b/src/lang/numeric/big_coefficient.h index df70905b7..ef14256e5 100644 --- a/src/lang/numeric/big_coefficient.h +++ b/src/lang/numeric/big_coefficient.h @@ -450,7 +450,7 @@ class BigCoefficient { trial_word * next_divisor_word > trial_remainder * BASE + normalized_dividend.words[position + divisor.length - 2]) { - trial_word--; + trial_word -= 1; trial_remainder += top_divisor_word; if (trial_remainder >= BASE) { break; @@ -473,7 +473,7 @@ class BigCoefficient { auto top_word = normalized_dividend.words[position + divisor.length]; if (top_word < borrow) { - trial_word--; + trial_word -= 1; std::uint64_t add_carry = 0; for (std::uint32_t index = 0; index < divisor.length; index++) { auto &word = normalized_dividend.words[position + index]; From b3aa02ab615d27feed65d1189b859c0ed2fd22b4 Mon Sep 17 00:00:00 2001 From: Juan Cruz Viotti Date: Fri, 7 Aug 2026 19:17:01 -0300 Subject: [PATCH 3/4] More Signed-off-by: Juan Cruz Viotti --- src/lang/numeric/big_coefficient.h | 59 ++++++++++++++++++++++++++++ src/lang/numeric/decimal.cc | 26 +++++++----- test/numeric/numeric_decimal_test.cc | 31 +++++++++++++++ 3 files changed, 106 insertions(+), 10 deletions(-) diff --git a/src/lang/numeric/big_coefficient.h b/src/lang/numeric/big_coefficient.h index ef14256e5..2c2a13ef2 100644 --- a/src/lang/numeric/big_coefficient.h +++ b/src/lang/numeric/big_coefficient.h @@ -513,6 +513,65 @@ class BigCoefficient { return {std::move(quotient), std::move(remainder)}; } + [[nodiscard]] auto multiply_modulo(const BigCoefficient &other, + const BigCoefficient &modulus) const + -> BigCoefficient { + auto product = this->multiply(other); + return product.divide_modulo(modulus).second; + } + + [[nodiscard]] static auto pow10_modulo(std::uint64_t power, + const BigCoefficient &modulus) + -> BigCoefficient { + BigCoefficient base{1}; + base.words[0] = 10; + base.length = 1; + base = base.divide_modulo(modulus).second; + + BigCoefficient result{1}; + result.words[0] = 1; + result.length = 1; + result = result.divide_modulo(modulus).second; + + auto remaining = power; + while (remaining > 0) { + if (remaining & 1) { + result = result.multiply_modulo(base, modulus); + } + + remaining >>= 1; + if (remaining > 0) { + base = base.multiply_modulo(base, modulus); + } + } + + return result; + } + + // The exponent difference between two decimal operands can reach billions, + // so scaling the dividend coefficient digit by digit before dividing would + // materialize gigabytes. Reducing the dividend first and folding the scale + // in through modular exponentiation keeps every intermediate bounded by the + // divisor size + [[nodiscard]] auto modulo_scaled(const BigCoefficient &divisor, + std::int64_t exponent_difference) const + -> BigCoefficient { + if (exponent_difference > 0) { + auto reduced = this->divide_modulo(divisor).second; + auto scale = pow10_modulo(static_cast(exponent_difference), + divisor); + return reduced.multiply_modulo(scale, divisor); + } + + if (exponent_difference < 0) { + auto scaled_divisor = divisor.multiply_pow10( + static_cast(-exponent_difference)); + return this->divide_modulo(scaled_divisor).second; + } + + return this->divide_modulo(divisor).second; + } + [[nodiscard]] static auto from_uint64(std::uint64_t value) -> BigCoefficient { if (value < BASE) { BigCoefficient result{1}; diff --git a/src/lang/numeric/decimal.cc b/src/lang/numeric/decimal.cc index 765f1192a..8c689aba0 100644 --- a/src/lang/numeric/decimal.cc +++ b/src/lang/numeric/decimal.cc @@ -1117,17 +1117,24 @@ auto Decimal::divisible_by(const Decimal &divisor) const -> bool { return static_cast(remaining % divisor_value) == 0; } + Decimal dividend_magnitude{*this}; + dividend_magnitude.flags_ = + static_cast(dividend_magnitude.flags_ & ~FLAG_SIGN); + Decimal divisor_magnitude{divisor}; + divisor_magnitude.flags_ = + static_cast(divisor_magnitude.flags_ & ~FLAG_SIGN); + if (dividend_magnitude < divisor_magnitude) { + return false; + } + auto dividend_big = coefficient_as_big(this->coefficient_, this->coefficient_high_, this->flags_); auto divisor_big = coefficient_as_big( divisor.coefficient_, divisor.coefficient_high_, divisor.flags_); - BigCoefficient::align_exponents(dividend_big, divisor_big, this->exponent_, - divisor.exponent_); - - auto [quotient, remainder] = dividend_big.divide_modulo(divisor_big); - - return remainder.is_zero(); + auto exponent_difference = + static_cast(this->exponent_) - divisor.exponent_; + return dividend_big.modulo_scaled(divisor_big, exponent_difference).is_zero(); } auto Decimal::same_quantum(const Decimal &other) const -> bool { @@ -2028,10 +2035,9 @@ auto Decimal::operator%=(const Decimal &other) -> Decimal & { auto divisor_big = coefficient_as_big(other.coefficient_, other.coefficient_high_, other.flags_); - BigCoefficient::align_exponents(dividend_big, divisor_big, this->exponent_, - other.exponent_); - - auto [quotient, remainder] = dividend_big.divide_modulo(divisor_big); + auto exponent_difference = + static_cast(this->exponent_) - other.exponent_; + auto remainder = dividend_big.modulo_scaled(divisor_big, exponent_difference); free_big_coefficient(this->coefficient_, this->flags_); store_big_result(this->coefficient_, this->coefficient_high_, this->flags_, diff --git a/test/numeric/numeric_decimal_test.cc b/test/numeric/numeric_decimal_test.cc index 4900cb866..8aeb3f70a 100644 --- a/test/numeric/numeric_decimal_test.cc +++ b/test/numeric/numeric_decimal_test.cc @@ -4448,6 +4448,37 @@ TEST(modulo_huge_divisor_returns_dividend) { EXPECT_EQ(dividend % divisor, sourcemeta::core::Decimal{"0.5"}); } +TEST(modulo_extreme_positive_dividend_exponent) { + const sourcemeta::core::Decimal dividend{"1e+2147483647"}; + const sourcemeta::core::Decimal divisor{1}; + EXPECT_TRUE((dividend % divisor).is_zero()); +} + +TEST(modulo_huge_exponent_multi_word_divisor) { + const sourcemeta::core::Decimal dividend{"1e+1000000"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_EQ((dividend % divisor).to_string(), "140735615018817496539760"); +} + +TEST(modulo_opposite_extreme_exponents) { + const sourcemeta::core::Decimal dividend{"1e+2147483647"}; + const sourcemeta::core::Decimal divisor{"1e-2147483648"}; + EXPECT_TRUE((dividend % divisor).is_zero()); +} + +TEST(divisible_by_huge_exponent_multi_word_divisor) { + const sourcemeta::core::Decimal dividend{"1e+1000000"}; + const sourcemeta::core::Decimal divisor{"1234567890123456789012345"}; + EXPECT_FALSE(dividend.divisible_by(divisor)); +} + +TEST(divisible_by_huge_exponent_power_of_ten_divisor) { + const sourcemeta::core::Decimal dividend{"1e+1000000"}; + const sourcemeta::core::Decimal divisor{ + "1000000000000000000000000000000000000000"}; + EXPECT_TRUE(dividend.divisible_by(divisor)); +} + TEST(add_exact_below_int64_minimum) { const sourcemeta::core::Decimal left{ std::numeric_limits::min()}; From 43bffc970635f4d595e0e13d8a9cf6c9cd1229a3 Mon Sep 17 00:00:00 2001 From: Juan Cruz Viotti Date: Fri, 7 Aug 2026 19:43:11 -0300 Subject: [PATCH 4/4] Fix Signed-off-by: Juan Cruz Viotti --- src/lang/numeric/big_coefficient.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/lang/numeric/big_coefficient.h b/src/lang/numeric/big_coefficient.h index 2c2a13ef2..2f147b83e 100644 --- a/src/lang/numeric/big_coefficient.h +++ b/src/lang/numeric/big_coefficient.h @@ -550,9 +550,12 @@ class BigCoefficient { // The exponent difference between two decimal operands can reach billions, // so scaling the dividend coefficient digit by digit before dividing would - // materialize gigabytes. Reducing the dividend first and folding the scale - // in through modular exponentiation keeps every intermediate bounded by the - // divisor size + // materialize gigabytes. When the dividend carries the larger exponent, + // reducing it first and folding the scale in through modular exponentiation + // keeps every intermediate bounded by the divisor size. When the divisor + // carries the larger exponent, its coefficient is scaled up by the full + // difference, so callers must first rule out dividends smaller in magnitude + // than the divisor, which bounds that scaling by the dividend digit count [[nodiscard]] auto modulo_scaled(const BigCoefficient &divisor, std::int64_t exponent_difference) const -> BigCoefficient {