Fix Decimal division hang and make additive arithmetic exact

src/lang/numeric/big_coefficient.h

@@ -393,54 +393,188 @@ 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<sourcemeta::core::uint128_t>(
-          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<std::uint64_t>(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<sourcemeta::core::uint128_t>(divisor.words[index]) *
+              normalizer +
+          normalize_carry;
+      normalized_divisor.words[index] =
+          static_cast<std::uint64_t>(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<sourcemeta::core::uint128_t>(this->words[index]) *
+              normalizer +
+          normalize_carry;
+      normalized_dividend.words[index] =
+          static_cast<std::uint64_t>(product % BASE);
+      normalize_carry = product / BASE;
+    }
+
+    normalized_dividend.words[this->length] =
+        static_cast<std::uint64_t>(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<sourcemeta::core::uint128_t>(
+              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 -= 1;
+        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<std::uint64_t>(subtrahend % BASE);
+        borrow = static_cast<std::uint64_t>(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 -= 1;
+        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<std::uint64_t>(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<std::uint64_t>(current / normalizer);
+      denormalize_carry = current % normalizer;
     }
 
     quotient.trim();
     remainder.trim();
     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. 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 {
+    if (exponent_difference > 0) {
+      auto reduced = this->divide_modulo(divisor).second;
+      auto scale = pow10_modulo(static_cast<std::uint64_t>(exponent_difference),
+                                divisor);
+      return reduced.multiply_modulo(scale, divisor);
+    }
+
+    if (exponent_difference < 0) {
+      auto scaled_divisor = divisor.multiply_pow10(
+          static_cast<std::uint32_t>(-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};

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<float>::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<double>::quiet_NaN();
+  }
+
   try {
     return std::stod(this->to_scientific_string());
   } catch (const std::out_of_range &) {
@@ -1103,17 +1117,24 @@ auto Decimal::divisible_by(const Decimal &divisor) const -> bool {
     return static_cast<std::uint64_t>(remaining % divisor_value) == 0;
   }
 
+  Decimal dividend_magnitude{*this};
+  dividend_magnitude.flags_ =
+      static_cast<std::uint8_t>(dividend_magnitude.flags_ & ~FLAG_SIGN);
+  Decimal divisor_magnitude{divisor};
+  divisor_magnitude.flags_ =
+      static_cast<std::uint8_t>(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<std::int64_t>(this->exponent_) - divisor.exponent_;
+  return dividend_big.modulo_scaled(divisor_big, exponent_difference).is_zero();
 }
 
 auto Decimal::same_quantum(const Decimal &other) const -> bool {
@@ -1791,8 +1812,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 +2012,37 @@ auto Decimal::operator%=(const Decimal &other) -> Decimal & {
     return *this;
   }
 
-  Decimal quotient{*this};
-  quotient /= other;
-
-  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<std::int32_t>(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<std::size_t>(number_of_digits - digits_to_remove));
-          auto old_sign =
-              static_cast<std::uint8_t>(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<std::uint8_t>(quotient.flags_ | old_sign);
-        }
+  // 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<std::uint8_t>(dividend_magnitude.flags_ & ~FLAG_SIGN);
+  Decimal divisor_magnitude{other};
+  divisor_magnitude.flags_ =
+      static_cast<std::uint8_t>(divisor_magnitude.flags_ & ~FLAG_SIGN);
+  if (dividend_magnitude < divisor_magnitude) {
+    return *this;
+  }
 
-      } else {
-        auto coefficient = quotient.coefficient_;
-        auto exponent = quotient.exponent_;
-        while (exponent < 0 && coefficient > 0) {
-          coefficient /= 10;
-          exponent++;
-        }
+  bool result_negative = (this->flags_ & FLAG_SIGN) != 0;
+  auto result_exponent = std::min(this->exponent_, other.exponent_);
 
-        if (exponent < 0) {
-          quotient = Decimal{};
-        } else {
-          quotient.coefficient_ = coefficient;
-          quotient.exponent_ = 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_);
 
-  Decimal product{quotient};
-  product *= other;
-  *this -= product;
+  auto exponent_difference =
+      static_cast<std::int64_t>(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_,
+                   std::move(remainder), result_negative);
+  this->exponent_ = result_exponent;
   return *this;
 }
 

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<std::int64_t>::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<std::int64_t>::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) {