Correctness fixes

lib/decimal.ex

@@ -984,7 +984,7 @@ defmodule Decimal do
 
         case integer_division(div_sign, coef1, exp1, coef2, exp2) do
           {:ok, result} ->
-            sub(num1, mult(num2, result))
+            exact_rem(num1, num2, result)
 
           {:error, error, reason, num} ->
             error(error, reason, num)
@@ -1065,7 +1065,7 @@ defmodule Decimal do
       true ->
         case integer_division(div_sign, coef1, exp1, coef2, exp2) do
           {:ok, result} ->
-            {result, sub(num1, mult(num2, result))}
+            {result, exact_rem(num1, num2, result)}
 
           {:error, error, reason, num} ->
             error(error, reason, {num, num})
@@ -1464,8 +1464,14 @@ defmodule Decimal do
       context(%Decimal{sign: 1, coef: coef, exp: exp}, [], false, ctx)
     else
       # otherwise the calculated root is inexact (but still meets precision),
-      # so use the root as `coef` and get the final exponent by shifting `exp`
-      context(%Decimal{sign: 1, coef: root, exp: exp - shift}, [], false, ctx)
+      # so use the root as `coef` and get the final exponent by shifting `exp`.
+      # The true root lies strictly beyond the truncated `root` on this branch
+      # (either the root itself is inexact, or an inexact down-shift dropped
+      # digits from an exact one), so the sticky bit is set: without it,
+      # rounding treats the guard digit as the entire discarded part, leaving
+      # directed modes short of the true root and turning "guard 5 with more
+      # beyond" into a false half-even tie.
+      context(%Decimal{sign: 1, coef: root, exp: exp - shift}, [], true, ctx)
     end
   end
 
@@ -2164,10 +2170,14 @@ defmodule Decimal do
     quo = Kernel.div(num, den)
     rem = num - quo * den
 
+    # Compare the doubled remainder so halving `den` cannot floor: with an
+    # exact division (`rem` 0) against `den` of 1, `den >>> 1` would be 0 and
+    # the exact quotient would fall through to the ties-to-even clause,
+    # rounding odd 53-bit significands away by one ULP.
     tmp =
-      case den >>> 1 do
-        den when rem > den -> quo + 1
-        den when rem < den -> quo
+      case rem <<< 1 do
+        rem2 when rem2 > den -> quo + 1
+        rem2 when rem2 < den -> quo
         _ when (quo &&& 1) === 1 -> quo + 1
         _ -> quo
       end
@@ -2431,6 +2441,26 @@ defmodule Decimal do
     }
   end
 
+  # The remainder is computed exactly: routing `num2 * quotient` through the
+  # public operations would round the product to the context precision before
+  # the subtraction, and a rounded product can cancel the true remainder
+  # entirely (for 34-digit operands the rounded product may equal `num1`).
+  # Only the final remainder goes through the context, like any result. The
+  # intermediates stay input-proportional: `integer_division/5` caps the
+  # quotient at `precision + 1` digits, which also bounds the exponent gap
+  # the alignment bridges.
+  #
+  # The quotient is floor(|num1| / |num2|), so the difference is non-negative
+  # and the remainder carries the dividend's sign - also when the remainder
+  # is zero, as IEEE 754 defines it.
+  defp exact_rem(%Decimal{} = num1, %Decimal{} = num2, %Decimal{coef: qcoef}) do
+    %Decimal{sign: sign1, coef: coef1, exp: exp1} = num1
+    %Decimal{coef: coef2, exp: exp2} = num2
+
+    {coef1, prod} = add_align(coef1, exp1, coef2 * qcoef, exp2)
+    context(%Decimal{sign: sign1, coef: coef1 - prod, exp: Kernel.min(exp1, exp2)})
+  end
+
   defp do_normalize(coef, exp) when coef >= @normalize_chunk_pow do
     case Kernel.rem(coef, @normalize_chunk_pow) do
       0 ->

test/decimal_test.exs

@@ -604,6 +604,34 @@ defmodule DecimalTest do
     end
   end
 
+  test "rem/2 and div_rem/2 compute the remainder exactly" do
+    # 34-digit operands whose divisor * quotient spans 67 digits: rounding
+    # that intermediate product to the context precision yields exactly the
+    # dividend, cancelling the true remainder of 3E-33 down to 0.
+    x = ~d"9999999999999999999999999999999999"
+    y = ~d"2.000000000000000000000000000000001"
+
+    assert Decimal.rem(x, y) == d(1, 3, -33)
+    assert Decimal.rem(~d"-9999999999999999999999999999999999", y) == d(-1, 3, -33)
+
+    {q, r} = Decimal.div_rem(x, y)
+    assert q == d(1, 4_999_999_999_999_999_999_999_999_999_999_997, 0)
+    assert r == d(1, 3, -33)
+
+    # the remainder is exact, so nothing may signal from the internals
+    flags = Context.get().flags
+    refute :inexact in flags
+    refute :rounded in flags
+
+    # a zero remainder takes the sign of the dividend, like any nonzero
+    # remainder does (and as IEEE 754 and Python's decimal define it)
+    assert Decimal.rem(~d"-4", ~d"2") == d(-1, 0, 0)
+
+    Context.with(%Context{precision: 5}, fn ->
+      assert Decimal.rem(~d"99999", ~d"2.0001") == d(1, 3, -4)
+    end)
+  end
+
   test "max/2" do
     assert Decimal.max(~d"0", ~d"0") == d(1, 0, 0)
     assert Decimal.max(~d"1", ~d"0") == d(1, 1, 0)
@@ -997,6 +1025,20 @@ defmodule DecimalTest do
     end)
   end
 
+  test "to_float/1 odd integers with 53-bit significands are exact" do
+    # 2^52 + 1, 2^53 - 1, and any odd integer between them are exactly
+    # representable as doubles and must convert without rounding.
+    assert Decimal.to_float(~d"4503599627370497") === 4_503_599_627_370_497.0
+    assert Decimal.to_float(~d"-4503599627370497") === -4_503_599_627_370_497.0
+    assert Decimal.to_float(~d"6004799503160661") === 6_004_799_503_160_661.0
+    assert Decimal.to_float(~d"9007199254740991") === 9_007_199_254_740_991.0
+
+    # Just above 2^53 the ULP is 2: odd values are genuine ties and must
+    # keep rounding to the even significand.
+    assert Decimal.to_float(~d"9007199254740993") === 9_007_199_254_740_992.0
+    assert Decimal.to_float(~d"9007199254740995") === 9_007_199_254_740_996.0
+  end
+
   test "round/3: special" do
     assert Decimal.round(~d"inf", 2, :down) == d(1, :inf, 0)
     assert Decimal.round(~d"nan", 2, :down) == d(1, :NaN, 0)
@@ -1120,6 +1162,40 @@ defmodule DecimalTest do
     end)
   end
 
+  test "sqrt/1 carries the discarded digits into rounding as a sticky bit" do
+    # An inexact root lies strictly beyond its truncated coefficient, so
+    # directed roundings must bump even when the guard digit is 0:
+    # sqrt(10) = 3.16227766016... > 3.16227766.
+    Context.with(%Context{precision: 9, rounding: :ceiling}, fn ->
+      assert Decimal.sqrt(~d"10") == d(1, 316_227_767, -8)
+      assert Decimal.sqrt(~d"2") == d(1, 141_421_357, -8)
+    end)
+
+    Context.with(%Context{precision: 9, rounding: :up}, fn ->
+      assert Decimal.sqrt(~d"10") == d(1, 316_227_767, -8)
+    end)
+
+    # :floor truncates a positive result regardless of discarded digits.
+    Context.with(%Context{precision: 9, rounding: :floor}, fn ->
+      assert Decimal.sqrt(~d"10") == d(1, 316_227_766, -8)
+    end)
+
+    # A guard digit of 5 with nonzero digits beyond it is not a tie:
+    # sqrt(1.57) = 1.25299... rounds up, not to the even neighbor.
+    Context.with(%Context{precision: 2, rounding: :half_even}, fn ->
+      assert Decimal.sqrt(~d"1.57") == d(1, 13, -1)
+    end)
+
+    # An inexact root signals :inexact, not just :rounded.
+    Context.with(%Context{precision: 9, rounding: :half_even}, fn ->
+      assert Decimal.sqrt(~d"10") == d(1, 316_227_766, -8)
+
+      flags = Context.get().flags
+      assert :inexact in flags
+      assert :rounded in flags
+    end)
+  end
+
   test "integer?/1" do
     assert Decimal.integer?(~d"1.0000")
     assert Decimal.integer?(~d"1")