From f1764cd97cce033fecaead94c5594e2612ad919f Mon Sep 17 00:00:00 2001 From: Rose Hogenson Date: Mon, 6 Jan 2025 04:40:07 -0800 Subject: Fix a rounding bug in int_to_float. My tiny brain doesn't really do well with this stuff. --- src/eval.rs | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) (limited to 'src/eval.rs') diff --git a/src/eval.rs b/src/eval.rs index d9fb184..beda22b 100644 --- a/src/eval.rs +++ b/src/eval.rs @@ -14,26 +14,31 @@ fn int_to_float(i: &BigInt) -> f64 { } let mut fraction; if exponent <= 52 { + // The number can be represented exactly without rounding. fraction = i.to_u64().unwrap() << 52 - exponent; } else { + // We need to round. let power = BigInt::from(1) << exponent - 52; - let (fraction_big, rem) = i.div_rem(&power); - fraction = fraction_big.to_u64().unwrap(); + let (rounded_fraction, rem) = i.div_rem(&power); + fraction = rounded_fraction.to_u64().unwrap(); let half = power >> 1; + // Round up? if rem > half || rem == half && fraction % 2 == 1 { fraction += 1; + // Check for overflow. if fraction == 1 << 53 { exponent += 1; + fraction >>= 1; if exponent > 1023 { return f64::INFINITY; } - fraction = 1 << 53; } } } - let u = (exponent + 1023) << 52 | fraction & !(1 << 52); - f64::from_bits(u) + // Build the float from exponent and fraction + // https://en.wikipedia.org/wiki/IEEE_754 + f64::from_bits(((exponent + 1023) << 52) + (fraction & !(1 << 52))) } #[derive(Debug, PartialEq)] @@ -301,6 +306,14 @@ mod tests { ); } + #[test] + fn int_to_float_round_up_overflow() { + assert_eq!( + int_to_float(&BigInt::from(0x1ff_ffff_ffff_fff9u64)), + 0x200_0000_0000_0000u64 as f64 + ); + } + #[test] fn int_to_float_round_down() { assert_eq!( -- cgit v1.3.1