Implement ties-to-even for Big[U]Int::to_{f32,f64}

BigUInt::to_{f32,f64} uses only the highest 64 bits as the mantissa,
discarding the rest, which gives an incorrect result when later
converted to float using IEEE-754's default nearest-ties-to-even
rounding.

Fixed to use correct round-to-odd for taking the mantissa, which means
setting the LSB to 1 if any of the lower (i.e. rounded away) bits are 1.

This also fixes the problem for BigInt since it really just uses BigUInt
under the hood.

src/biguint/convert.rs

@@ -287,19 +287,31 @@ fn high_bits_to_u64(v: &BigUint) -> u64 {
                 let digit_bits = (bits - 1) % u64::from(big_digit::BITS) + 1;
                 let bits_want = Ord::min(64 - ret_bits, digit_bits);
 
-                if bits_want != 64 {
-                    ret <<= bits_want;
+                if bits_want != 0 {
+                    if bits_want != 64 {
+                        ret <<= bits_want;
+                    }
+                    // XXX Conversion is useless if already 64-bit.
+                    #[allow(clippy::useless_conversion)]
+                    let d0 = u64::from(*d) >> (digit_bits - bits_want);
+                    ret |= d0;
                 }
-                // XXX Conversion is useless if already 64-bit.
-                #[allow(clippy::useless_conversion)]
-                let d0 = u64::from(*d) >> (digit_bits - bits_want);
-                ret |= d0;
-                ret_bits += bits_want;
-                bits -= bits_want;
 
-                if ret_bits == 64 {
-                    break;
+                // Implement round-to-odd: If any lower bits are 1, set LSB to 1
+                // so that rounding again to floating point value using
+                // nearest-ties-to-even is correct.
+                //
+                // See: https://en.wikipedia.org/wiki/Rounding#Rounding_to_prepare_for_shorter_precision
+
+                if digit_bits - bits_want != 0 {
+                    // XXX Conversion is useless if already 64-bit.
+                    #[allow(clippy::useless_conversion)]
+                    let masked = u64::from(*d) << (64 - (digit_bits - bits_want) as u32);
+                    ret |= (masked != 0) as u64;
                 }
+
+                ret_bits += bits_want;
+                bits -= bits_want;
             }
 
             ret

tests/bigint.rs

@@ -421,6 +421,13 @@ fn test_convert_f32() {
         b <<= 1;
     }
 
+    // test correct ties-to-even rounding
+    let weird: i128 = (1i128 << 100) + (1i128 << (100 - f32::MANTISSA_DIGITS));
+    assert_ne!(weird as f32, (weird + 1) as f32);
+
+    assert_eq!(BigInt::from(weird).to_f32(), Some(weird as f32));
+    assert_eq!(BigInt::from(weird + 1).to_f32(), Some((weird + 1) as f32));
+
     // rounding
     assert_eq!(
         BigInt::from_f32(-f32::consts::PI),
@@ -511,6 +518,13 @@ fn test_convert_f64() {
         b <<= 1;
     }
 
+    // test correct ties-to-even rounding
+    let weird: i128 = (1i128 << 100) + (1i128 << (100 - f64::MANTISSA_DIGITS));
+    assert_ne!(weird as f64, (weird + 1) as f64);
+
+    assert_eq!(BigInt::from(weird).to_f64(), Some(weird as f64));
+    assert_eq!(BigInt::from(weird + 1).to_f64(), Some((weird + 1) as f64));
+
     // rounding
     assert_eq!(
         BigInt::from_f64(-f64::consts::PI),

tests/biguint.rs

@@ -646,6 +646,13 @@ fn test_convert_f32() {
         b <<= 1;
     }
 
+    // test correct ties-to-even rounding
+    let weird: i128 = (1i128 << 100) + (1i128 << (100 - f32::MANTISSA_DIGITS));
+    assert_ne!(weird as f32, (weird + 1) as f32);
+
+    assert_eq!(BigInt::from(weird).to_f32(), Some(weird as f32));
+    assert_eq!(BigInt::from(weird + 1).to_f32(), Some((weird + 1) as f32));
+
     // rounding
     assert_eq!(BigUint::from_f32(-1.0), None);
     assert_eq!(BigUint::from_f32(-0.99999), Some(BigUint::zero()));
@@ -722,6 +729,13 @@ fn test_convert_f64() {
         b <<= 1;
     }
 
+    // test correct ties-to-even rounding
+    let weird: i128 = (1i128 << 100) + (1i128 << (100 - f64::MANTISSA_DIGITS));
+    assert_ne!(weird as f64, (weird + 1) as f64);
+
+    assert_eq!(BigInt::from(weird).to_f64(), Some(weird as f64));
+    assert_eq!(BigInt::from(weird + 1).to_f64(), Some((weird + 1) as f64));
+
     // rounding
     assert_eq!(BigUint::from_f64(-1.0), None);
     assert_eq!(BigUint::from_f64(-0.99999), Some(BigUint::zero()));