diff --git a/float/CHANGELOG.md b/float/CHANGELOG.md index f3c99017..4d450cc2 100644 --- a/float/CHANGELOG.md +++ b/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 diff --git a/float/src/convert.rs b/float/src/convert.rs index 671f573a..a73d1e1e 100644 --- a/float/src/convert.rs +++ b/float/src/convert.rs @@ -433,15 +433,7 @@ impl FBig { /// ``` #[inline] pub fn to_f32(&self) -> Rounded { - if self.repr.is_infinite() { - return Inexact(self.sign() * f32::INFINITY, Rounding::NoOp); - } - - let context = Context::::new(24); - context - .convert_base::(self.repr.clone(), None) - .and_then(|v| context.repr_round_ref(&v)) - .and_then(|v| v.into_f32_internal()) + Context::::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 FBig { /// ``` #[inline] pub fn to_f64(&self) -> Rounded { - if self.repr.is_infinite() { - return Inexact(self.sign() * f64::INFINITY, Rounding::NoOp); - } - - let context = Context::::new(53); - context - .convert_base::(self.repr.clone(), None) - .and_then(|v| context.repr_round_ref(&v)) - .and_then(|v| v.into_f64_internal()) + Context::::convert_to_f64(self.repr.clone()) } } @@ -488,7 +472,67 @@ fn converted_overflow_repr(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(repr: Repr, width: usize) -> Repr<2> { + match Context::::new(width).convert_base::(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 Context { + // 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(repr: Repr) -> Rounded { + if repr.is_infinite() { + return Inexact(repr.sign() * f64::INFINITY, Rounding::NoOp); + } + let odd = convert_base_odd::(repr, 60); + let bits = significand_bits(&odd, 53, -1074); + Context::::new(bits) + .repr_round(odd) + .and_then(|v| v.into_f64_internal()) + } + + // [convert_to_f64] for f32. + fn convert_to_f32(repr: Repr) -> Rounded { + if repr.is_infinite() { + return Inexact(repr.sign() * f32::INFINITY, Rounding::NoOp); + } + let odd = convert_base_odd::(repr, 32); + let bits = significand_bits(&odd, 24, -149); + Context::::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( @@ -588,9 +632,22 @@ impl Context { let signif = repr.significand * Repr::::BASE.pow(repr.exponent as usize); Exact(Repr::new(signif, 0)) } else { - let num: Repr = Repr::new(repr.significand, 0); let den: Repr = Repr::new(Repr::::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 = 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| Repr { significand: r.significand, @@ -680,15 +737,7 @@ impl Repr { /// ``` #[inline] pub fn to_f32(&self) -> Rounded { - if self.is_infinite() { - return Inexact(self.sign() * f32::INFINITY, Rounding::NoOp); - } - - let context = Context::::new(24); - context - .convert_base::(self.clone(), None) - .and_then(|v| context.repr_round_ref(&v)) - .and_then(|v| v.into_f32_internal()) + Context::::convert_to_f32(self.clone()) } // this method requires that the representation is already rounded to 53 binary bits @@ -736,15 +785,7 @@ impl Repr { /// ``` #[inline] pub fn to_f64(&self) -> Rounded { - if self.is_infinite() { - return Inexact(self.sign() * f64::INFINITY, Rounding::NoOp); - } - - let context = Context::::new(53); - context - .convert_base::(self.clone(), None) - .and_then(|v| context.repr_round_ref(&v)) - .and_then(|v| v.into_f64_internal()) + Context::::convert_to_f64(self.clone()) } /// Convert the float number representation to a [IBig]. diff --git a/float/tests/convert.rs b/float/tests/convert.rs index 74ccfc66..89a8ea20 100644 --- a/float/tests/convert.rs +++ b/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 + δ + 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::::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 + δ + 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::::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}" + ); + } +}