Fix incorrect rounding with subnormal/zero results of float multiplication

Author: beetrees

src/float/mul.rs

@@ -149,18 +149,13 @@ where
         }
 
         // Otherwise, shift the significand of the result so that the round
-        // bit is the high bit of productLo.
-        if shift < bits {
-            let sticky = product_low << (bits - shift);
-            product_low = product_high << (bits - shift) | product_low >> shift | sticky;
-            product_high >>= shift;
-        } else if shift < (2 * bits) {
-            let sticky = product_high << (2 * bits - shift) | product_low;
-            product_low = product_high >> (shift - bits) | sticky;
-            product_high = zero;
-        } else {
-            product_high = zero;
-        }
+        // bit is the high bit of `product_low`.
+        // Ensure one of the non-highest bits in `product_low` is set if the shifted out bit are
+        // not all zero so that the result is correctly rounded below.
+        let sticky = product_low << (bits - shift) != zero;
+        product_low =
+            product_high << (bits - shift) | product_low >> shift | (sticky as u32).cast();
+        product_high >>= shift;
     } else {
         // Result is normal before rounding; insert the exponent.
         product_high &= significand_mask;

testcrate/tests/mul.rs

@@ -107,14 +107,11 @@ macro_rules! float_mul {
                 fuzz_float_2(N, |x: $f, y: $f| {
                     let mul0 = apfloat_fallback!($f, $apfloat_ty, $sys_available, Mul::mul, x, y);
                     let mul1: $f = $fn(x, y);
-                    // multiplication of subnormals is not currently handled
-                    if !(Float::is_subnormal(mul0) || Float::is_subnormal(mul1)) {
-                        if !Float::eq_repr(mul0, mul1) {
-                            panic!(
-                                "{}({:?}, {:?}): std: {:?}, builtins: {:?}",
-                                stringify!($fn), x, y, mul0, mul1
-                            );
-                        }
+                    if !Float::eq_repr(mul0, mul1) {
+                        panic!(
+                            "{}({:?}, {:?}): std: {:?}, builtins: {:?}",
+                            stringify!($fn), x, y, mul0, mul1
+                        );
                     }
                 });
             }