Round float at the decimal level to match CPython's _Py_dg_dtoa (#7761)

* Round float at the decimal level to match CPython's _Py_dg_dtoa

CPython's `float.__round__` (Objects/floatobject.c) routes through
`_Py_dg_dtoa` and rounds at the decimal level. The previous
`round_float_digits` multiplied by 10**ndigits and rounded at the
IEEE 754 binary level, which diverges for values that aren't exactly
representable. For example, 2.675 stores as 2.67499...; dtoa correctly
rounds it down to 2.67, but `(2.675 * 100.0).round() / 100.0` lands on
2.68 because the multiplication produces a phantom 267.5 tie that
round-half-to-even snaps up.

Rust's `{:.*}` float formatting uses dtoa-style algorithms (Grisu3 +
Dragon4 fallback) and matches CPython's `_Py_dg_dtoa` byte-for-byte.
Replace the multiply-then-round step with `format!` + `parse` for
ndigits >= 0. The ndigits < 0 path is unchanged because dividing
typical inputs by 10**|ndigits| produces genuine ties rather than
synthesizing them.

Verified byte-identical with CPython 3.14.4 over a 108-case random
fuzz plus targeted half-tie probes. Unmasks
`test_float.RoundTestCase.test_matches_float_format` and
`test_previous_round_bugs`.

* Use #[expect] with reason for float_cmp suppression

Co-authored-by: ShaharNaveh <50263213+ShaharNaveh@users.noreply.github.com>

---------

Co-authored-by: ShaharNaveh <50263213+ShaharNaveh@users.noreply.github.com>

Lib/test/test_float.py

@@ -965,7 +965,6 @@ class RoundTestCase(unittest.TestCase, FloatsAreIdenticalMixin):
         self.assertRaises(OverflowError, round, 1.6e308, -308)
         self.assertRaises(OverflowError, round, -1.7e308, -308)
 
-    @unittest.expectedFailure  # TODO: RUSTPYTHON; AssertionError: 56294995342131.51 != 56294995342131.5
     @unittest.skipUnless(getattr(sys, 'float_repr_style', '') == 'short',
                          "applies only when using short float repr style")
     def test_previous_round_bugs(self):
@@ -984,7 +983,6 @@ class RoundTestCase(unittest.TestCase, FloatsAreIdenticalMixin):
         self.assertEqual(round(85.0, -1), 80.0)
         self.assertEqual(round(95.0, -1), 100.0)
 
-    @unittest.expectedFailure  # TODO: RUSTPYTHON; AssertionError: 0.01 != 0.0
     @unittest.skipUnless(getattr(sys, 'float_repr_style', '') == 'short',
                          "applies only when using short float repr style")
     def test_matches_float_format(self):

crates/common/src/float_ops.rs

@@ -1,6 +1,6 @@
 use core::f64;
 use malachite_bigint::{BigInt, ToBigInt};
-use num_traits::{Float, Signed, ToPrimitive, Zero};
+use num_traits::{Signed, ToPrimitive};
 
 pub const fn decompose_float(value: f64) -> (f64, i32) {
     if 0.0 == value {
@@ -208,64 +208,53 @@ pub fn ulp(x: f64) -> f64 {
 }
 
 pub fn round_float_digits(x: f64, ndigits: i32) -> Option<f64> {
-    let float = if ndigits.is_zero() {
-        let fract = x.fract();
-        if (fract.abs() - 0.5).abs() < f64::EPSILON {
-            if x.trunc() % 2.0 == 0.0 {
-                x - fract
-            } else {
-                x + fract
-            }
-        } else {
-            x.round()
-        }
+    // Mirror CPython's `float.__round__` (Objects/floatobject.c), which uses
+    // `_Py_dg_dtoa` to round at the decimal level. Multiplying by 10**ndigits
+    // and rounding at the IEEE 754 binary level diverges for values that
+    // aren't exactly representable: 2.675 stores as 2.67499..., which dtoa
+    // correctly rounds down to 2.67, but `(2.675 * 100.0).round() / 100.0`
+    // lands on 2.68 because the multiplication produces a phantom 267.5 tie.
+    // Rust's `{:.*}` float formatting uses dtoa-style algorithms and matches
+    // CPython's `_Py_dg_dtoa` byte-for-byte.
+    if !x.is_finite() {
+        return Some(x);
+    }
+
+    const NDIGITS_MAX: i32 =
+        ((f64::MANTISSA_DIGITS as i32 - f64::MIN_EXP) as f64 * f64::consts::LOG10_2) as i32;
+    const NDIGITS_MIN: i32 = -(((f64::MAX_EXP + 1) as f64 * f64::consts::LOG10_2) as i32);
+
+    if ndigits > NDIGITS_MAX {
+        return Some(x);
+    }
+    if ndigits < NDIGITS_MIN {
+        return Some(0.0f64.copysign(x));
+    }
+
+    let result: f64 = if ndigits >= 0 {
+        let s = format!("{:.*}", ndigits as usize, x);
+        s.parse().ok()?
     } else {
-        const NDIGITS_MAX: i32 =
-            ((f64::MANTISSA_DIGITS as i32 - f64::MIN_EXP) as f64 * f64::consts::LOG10_2) as i32;
-        const NDIGITS_MIN: i32 = -(((f64::MAX_EXP + 1) as f64 * f64::consts::LOG10_2) as i32);
-        if ndigits > NDIGITS_MAX {
-            x
-        } else if ndigits < NDIGITS_MIN {
-            0.0f64.copysign(x)
+        // ndigits < 0: divide-then-round avoids the phantom-tie problem
+        // because dividing typical inputs by 10**|ndigits| produces genuine
+        // half-integer ties rather than synthesizing them.
+        let pow1 = 10.0f64.powi(-ndigits);
+        let y = x / pow1;
+        let z = y.round();
+        #[expect(
+            clippy::float_cmp,
+            reason = "exact half-tie detection for banker's rounding correction"
+        )]
+        let z = if (y - z).abs() == 0.5 {
+            2.0 * (y / 2.0).round()
         } else {
-            let (y, pow1, pow2) = if ndigits >= 0 {
-                // according to cpython: pow1 and pow2 are each safe from overflow, but
-                //                       pow1*pow2 ~= pow(10.0, ndigits) might overflow
-                let (pow1, pow2) = if ndigits > 22 {
-                    (10.0.powf((ndigits - 22) as f64), 1e22)
-                } else {
-                    (10.0.powf(ndigits as f64), 1.0)
-                };
-                let y = (x * pow1) * pow2;
-                if !y.is_finite() {
-                    return Some(x);
-                }
-                (y, pow1, Some(pow2))
-            } else {
-                let pow1 = 10.0.powf((-ndigits) as f64);
-                (x / pow1, pow1, None)
-            };
-            let z = y.round();
-            #[allow(clippy::float_cmp)]
-            let z = if (y - z).abs() == 0.5 {
-                2.0 * (y / 2.0).round()
-            } else {
-                z
-            };
-            let z = if let Some(pow2) = pow2 {
-                // ndigits >= 0
-                (z / pow2) / pow1
-            } else {
-                z * pow1
-            };
-
-            if !z.is_finite() {
-                // overflow
-                return None;
-            }
-
             z
-        }
+        };
+        z * pow1
     };
-    Some(float)
+
+    if !result.is_finite() {
+        return None;
+    }
+    Some(result)
 }

extra_tests/snippets/builtin_round.py

@@ -1,3 +1,5 @@
+import math
+
 from testutils import assert_raises
 
 assert round(1.2) == 1
@@ -43,3 +45,51 @@ class X:
 assert round(X(), 1) == 1.1
 assert round(X(), None) == 1.1
 assert round(X()) == 1.1
+
+
+# Banker's rounding at the decimal level: values like 2.675 store as
+# 2.67499... in IEEE 754, so a multiply-then-round implementation creates
+# a phantom 267.5 tie that round-half-to-even snaps up to 2.68. CPython's
+# `_Py_dg_dtoa` rounds at the decimal level and returns 2.67.
+assert round(2.675, 2) == 2.67
+assert round(2.685, 2) == 2.69
+assert round(-2.675, 2) == -2.67
+assert round(0.05, 1) == 0.1
+assert round(0.15, 1) == 0.1
+assert round(0.35, 1) == 0.3
+assert round(0.45, 1) == 0.5
+assert round(0.65, 1) == 0.7
+assert round(0.85, 1) == 0.8
+assert round(0.95, 1) == 0.9
+assert round(0.645, 2) == 0.65
+assert round(0.665, 2) == 0.67
+assert round(0.685, 2) == 0.69
+assert round(0.695, 2) == 0.69
+assert round(1.685, 2) == 1.69
+assert round(3.745, 2) == 3.75
+
+# Exact-halfway ties at integer level use banker's (round-half-to-even).
+assert round(0.5, 0) == 0.0
+assert round(1.5, 0) == 2.0
+assert round(2.5, 0) == 2.0
+
+# Negative ndigits uses divide-then-round; banker's still applies.
+assert round(1235, -1) == 1240
+assert round(1245, -1) == 1240
+assert round(1234.5, -1) == 1230.0
+assert round(150, -2) == 200
+assert round(250, -2) == 200
+
+# NaN and infinities round to themselves with no error when ndigits is given.
+assert math.isnan(round(float("nan"), 2))
+assert round(float("inf"), 2) == float("inf")
+assert round(float("-inf"), 2) == float("-inf")
+
+# Signed zero is preserved through rounding.
+assert math.copysign(1.0, round(-0.0, 2)) == -1.0
+assert math.copysign(1.0, round(0.0, 2)) == 1.0
+
+# Out-of-range ndigits short-circuits without overflow.
+assert round(1.0, 1000) == 1.0
+assert round(1.0, -1000) == 0.0
+assert round(1.7976931348623157e308, 0) == 1.7976931348623157e308