Round to_f64/to_f32 once, to the subnormal-aware precision

float/CHANGELOG.md

@@ -7,6 +7,16 @@
   significand (the `const` counterpart of `Repr::new`). `FBig::from_parts_const` now delegates to
   it, and the complex literal macro uses it.
 
+### Fix
+- `to_f64`/`to_f32` now round the source once, directly to the target's precision at its own
+  magnitude (fewer than 53/24 bits for subnormals), instead of through a fixed 53/24-bit
+  intermediate that re-rounds into the subnormal grid. This removes a 1-ULP double-rounding error
+  on subnormal values that sit just past a subnormal halfway.
+- `to_f64`/`to_f32` no longer panic in debug builds (nor silently double-round in release) on
+  high-precision inputs: the base-changing conversion now pre-shrinks the source significand before
+  dividing, upholding `repr_div`'s dividend-width contract (as `Context::div` already does) instead
+  of feeding an oversized dividend into the division.
+
 ## 0.5.0
 
 ### Add

float/src/convert.rs

@@ -433,15 +433,7 @@ impl<R: Round, const B: Word> FBig<R, B> {
     /// ```
     #[inline]
     pub fn to_f32(&self) -> Rounded<f32> {
-        if self.repr.is_infinite() {
-            return Inexact(self.sign() * f32::INFINITY, Rounding::NoOp);
-        }
-
-        let context = Context::<R>::new(24);
-        context
-            .convert_base::<B, 2>(self.repr.clone(), None)
-            .and_then(|v| context.repr_round_ref(&v))
-            .and_then(|v| v.into_f32_internal())
+        Context::<R>::convert_to_f32(self.repr.clone())
     }
 
     /// Convert the float number to [f64] with the rounding mode associated with the type.
@@ -460,15 +452,7 @@ impl<R: Round, const B: Word> FBig<R, B> {
     /// ```
     #[inline]
     pub fn to_f64(&self) -> Rounded<f64> {
-        if self.repr.is_infinite() {
-            return Inexact(self.sign() * f64::INFINITY, Rounding::NoOp);
-        }
-
-        let context = Context::<R>::new(53);
-        context
-            .convert_base::<B, 2>(self.repr.clone(), None)
-            .and_then(|v| context.repr_round_ref(&v))
-            .and_then(|v| v.into_f64_internal())
+        Context::<R>::convert_to_f64(self.repr.clone())
     }
 }
 
@@ -488,7 +472,67 @@ fn converted_overflow_repr<const NewB: Word>(large: bool, sign: Sign) -> Rounded
     )
 }
 
+/// Number of significant bits a binary float format keeps for a value whose most-significant bit
+/// sits at position `msb`: `max_bits` across the normal range, but fewer for subnormals, whose
+/// spacing is fixed at `2^subnormal_exp` (e.g. `2^-1074` for f64, `2^-149` for f32). Rounding the
+/// source straight to this width lets the bit-encoding step avoid a second rounding, which would
+/// otherwise double-round subnormals.
+fn significand_bits(v: &Repr<2>, max_bits: usize, subnormal_exp: isize) -> usize {
+    if v.significand.is_zero() {
+        return max_bits;
+    }
+    let msb = v.exponent + v.digits() as isize - 1;
+    (msb - subnormal_exp + 1).clamp(1, max_bits as isize) as usize
+}
+
+/// Convert `repr` to base 2 and truncate to `width` significant bits, forcing the lowest kept bit
+/// to 1 whenever the tail is nonzero (round-to-odd). Rounding this down to any width up to
+/// `width - 2` reproduces the correctly-rounded value for every rounding mode, so the two-step
+/// "convert, then round to the final width" cannot double-round. `width` is fixed and generous, so
+/// the base-conversion logarithm stays accurate even when the final width is tiny (deep subnormals).
+#[allow(non_upper_case_globals)]
+fn convert_base_odd<const B: Word>(repr: Repr<B>, width: usize) -> Repr<2> {
+    match Context::<Zero>::new(width).convert_base::<B, 2>(repr, None) {
+        Exact(v) => v,
+        Inexact(v, _) if v.significand.is_zero() => v,
+        Inexact(v, _) => {
+            let shift = width - v.digits();
+            let (sign, mut mag) = v.significand.into_parts();
+            mag <<= shift;
+            mag.set_bit(0);
+            Repr::new(IBig::from_parts(sign, mag), v.exponent - shift as isize)
+        }
+    }
+}
+
 impl<R: Round> Context<R> {
+    // Convert `repr` (base B) to the nearest f64 under this context's rounding mode. A generous
+    // round-to-odd base conversion is rounded once to the target's precision at its own magnitude
+    // (fewer than 53 bits for subnormals), so `into_f64_internal` re-rounds nothing — which would
+    // otherwise double-round subnormals. Handles a source significand of any size.
+    fn convert_to_f64<const B: Word>(repr: Repr<B>) -> Rounded<f64> {
+        if repr.is_infinite() {
+            return Inexact(repr.sign() * f64::INFINITY, Rounding::NoOp);
+        }
+        let odd = convert_base_odd::<B>(repr, 60);
+        let bits = significand_bits(&odd, 53, -1074);
+        Context::<R>::new(bits)
+            .repr_round(odd)
+            .and_then(|v| v.into_f64_internal())
+    }
+
+    // [convert_to_f64] for f32.
+    fn convert_to_f32<const B: Word>(repr: Repr<B>) -> Rounded<f32> {
+        if repr.is_infinite() {
+            return Inexact(repr.sign() * f32::INFINITY, Rounding::NoOp);
+        }
+        let odd = convert_base_odd::<B>(repr, 32);
+        let bits = significand_bits(&odd, 24, -149);
+        Context::<R>::new(bits)
+            .repr_round(odd)
+            .and_then(|v| v.into_f32_internal())
+    }
+
     // Convert the [Repr] from base B to base NewB, with the precision under the target base from this context.
     #[allow(non_upper_case_globals)]
     fn convert_base<const B: Word, const NewB: Word>(
@@ -588,9 +632,22 @@ impl<R: Round> Context<R> {
                 let signif = repr.significand * Repr::<B>::BASE.pow(repr.exponent as usize);
                 Exact(Repr::new(signif, 0))
             } else {
-                let num: Repr<NewB> = Repr::new(repr.significand, 0);
                 let den: Repr<NewB> =
                     Repr::new(Repr::<B>::BASE.pow(-repr.exponent as usize).into(), 0);
+                // repr_div requires the dividend to be no wider than `precision + divisor`, so
+                // pre-shrink the significand the same way Context::div does — the caller, not
+                // the kernel, is responsible for bounding the dividend. Rounding it to
+                // `den.digits() + precision` preserves enough information for the division to
+                // be correctly rounded at `precision`.
+                let num: Repr<NewB> = Repr::new(repr.significand, 0);
+                let num =
+                    if !num.is_pos_zero() && num.digits_ub() > den.digits_lb() + self.precision {
+                        Self::new(den.digits() + self.precision)
+                            .repr_round_ref(&num)
+                            .value()
+                    } else {
+                        num
+                    };
                 match self.repr_div(num, den) {
                     Ok(v) => v.map(|r: Repr<NewB>| Repr {
                         significand: r.significand,
@@ -680,15 +737,7 @@ impl<const B: Word> Repr<B> {
     /// ```
     #[inline]
     pub fn to_f32(&self) -> Rounded<f32> {
-        if self.is_infinite() {
-            return Inexact(self.sign() * f32::INFINITY, Rounding::NoOp);
-        }
-
-        let context = Context::<HalfEven>::new(24);
-        context
-            .convert_base::<B, 2>(self.clone(), None)
-            .and_then(|v| context.repr_round_ref(&v))
-            .and_then(|v| v.into_f32_internal())
+        Context::<HalfEven>::convert_to_f32(self.clone())
     }
 
     // this method requires that the representation is already rounded to 53 binary bits
@@ -736,15 +785,7 @@ impl<const B: Word> Repr<B> {
     /// ```
     #[inline]
     pub fn to_f64(&self) -> Rounded<f64> {
-        if self.is_infinite() {
-            return Inexact(self.sign() * f64::INFINITY, Rounding::NoOp);
-        }
-
-        let context = Context::<HalfEven>::new(53);
-        context
-            .convert_base::<B, 2>(self.clone(), None)
-            .and_then(|v| context.repr_round_ref(&v))
-            .and_then(|v| v.into_f64_internal())
+        Context::<HalfEven>::convert_to_f64(self.clone())
     }
 
     /// Convert the float number representation to a [IBig].

float/tests/convert.rs

@@ -3,11 +3,12 @@ use core::str::FromStr;
 use dashu_base::{Approximation::*, ConversionError::*};
 use dashu_float::{
     round::{
-        mode::{HalfAway, Zero},
+        mode::{HalfAway, HalfEven, Zero},
         Rounding::*,
     },
     DBig, FBig,
 };
+use dashu_int::{IBig, UBig};
 
 mod helper_macros;
 
@@ -657,3 +658,102 @@ fn test_dbig_to_f32() {
     assert_eq!(DBig::INFINITY.to_f32(), Inexact(f32::INFINITY, NoOp));
     assert_eq!(DBig::NEG_INFINITY.to_f32(), Inexact(f32::NEG_INFINITY, NoOp));
 }
+
+// A subnormal that sits just past a subnormal halfway must round straight to the
+// target's exponent-dependent precision. Rounding through a fixed 53-bit intermediate
+// drops the excess that lifts the value above the halfway, lands on the halfway, and
+// the final encoding rounds to even -- one ULP below the correct result. The exact
+// value (2m+1)/2^1075 +/- 1/2^(1075+j) is such a halfway nudged by a bit far below the
+// 53rd; correctly rounded it is m+1 (nudged up) or m (nudged down), and the subnormal
+// k*2^-1074 has bit pattern exactly k. Built in both base 2 and base 10 so the
+// base-changing conversion path is exercised too.
+#[test]
+fn test_to_f64_subnormal_halfway() {
+    fn check(m: u64, j: u32) {
+        let scale = (1075 + j) as usize;
+        let core = IBig::from(2 * m + 1) << j as usize;
+        let five: IBig = UBig::from(5u8).pow(scale).into();
+        for (delta, want) in [(IBig::ONE, m + 1), (-IBig::ONE, m)] {
+            // magnitude, always positive; cover both signs (negation flips only the sign bit)
+            let mag = &core + &delta;
+            for negative in [false, true] {
+                let sig = if negative { -mag.clone() } else { mag.clone() };
+                let want = if negative { want | (1u64 << 63) } else { want };
+                let base2 = FBig::<HalfEven, 2>::from_parts(sig.clone(), -(scale as isize));
+                assert_eq!(base2.to_f64().value().to_bits(), want, "base 2, m={m}, neg={negative}");
+                let base10 = DBig::from_parts(&sig * &five, -(scale as isize));
+                assert_eq!(
+                    base10.to_f64().value().to_bits(),
+                    want,
+                    "base 10, m={m}, neg={negative}"
+                );
+            }
+        }
+    }
+    // m spans the subnormal significand width; j puts the nudge well below the 53rd bit
+    // (but within the conversion's working precision), where a second rounding would
+    // otherwise collapse the value onto the halfway.
+    for m in [
+        3,
+        21,
+        (1 << 8) + 5,
+        (1 << 20) + 3,
+        (1 << 33) + 7,
+        (1 << 40) + 9,
+    ] {
+        check(m, 60);
+    }
+}
+
+// [test_to_f64_subnormal_halfway] for f32 (subnormal spacing 2^-149, 24-bit mantissa).
+#[test]
+fn test_to_f32_subnormal_halfway() {
+    fn check(m: u32, j: u32) {
+        let scale = (150 + j) as usize;
+        let core = IBig::from(2 * m + 1) << j as usize;
+        let five: IBig = UBig::from(5u8).pow(scale).into();
+        for (delta, want) in [(IBig::ONE, m + 1), (-IBig::ONE, m)] {
+            // magnitude, always positive; cover both signs (negation flips only the sign bit)
+            let mag = &core + &delta;
+            for negative in [false, true] {
+                let sig = if negative { -mag.clone() } else { mag.clone() };
+                let want = if negative { want | (1u32 << 31) } else { want };
+                let base2 = FBig::<HalfEven, 2>::from_parts(sig.clone(), -(scale as isize));
+                assert_eq!(base2.to_f32().value().to_bits(), want, "base 2, m={m}, neg={negative}");
+                let base10 = DBig::from_parts(&sig * &five, -(scale as isize));
+                assert_eq!(
+                    base10.to_f32().value().to_bits(),
+                    want,
+                    "base 10, m={m}, neg={negative}"
+                );
+            }
+        }
+    }
+    for m in [3, 21, (1 << 10) + 5, (1 << 18) + 7] {
+        check(m, 30);
+    }
+}
+
+// High-precision decimals feed an oversized significand into the base-changing
+// division; this tripped a debug assertion (and double-rounded in release). Each must
+// now convert without panicking to the correctly rounded f64 (oracle: Python float()).
+#[test]
+fn test_to_f64_high_precision() {
+    let cases = [
+        ("123456789012345678.9012345678901", 0x437b69b4ba630f35u64),
+        ("3915263378237002511617337316730e-19", 0x4256ca327347ecd1),
+        ("1234567890123456789012345678901e-5", 0x45246c993044fd55),
+        ("9999999999999999999999999999999e-3", 0x45c027e72f1f1281),
+        ("27182818284590452353602874713526e-13", 0x43c2dca375e059b1),
+    ];
+    for (s, want) in cases {
+        assert_eq!(DBig::from_str(s).unwrap().to_f64().value().to_bits(), want, "{s}");
+        // the negation flips only the sign bit; cover the sign through the over-wide-quotient path
+        let neg = format!("-{s}");
+        assert_eq!(
+            DBig::from_str(&neg).unwrap().to_f64().value().to_bits(),
+            want | (1u64 << 63),
+            "{neg}"
+        );
+    }
+}