Improve precision of Duration-float operations

Rather than convert Duration to a float in order to multiply or divide
by a floating point number, which entails a loss of precision, we
instead represent both operands as u128 nanoseconds and perform integer
arithmetic thereupon.

Author: eggyal

library/core/src/time.rs

@@ -1020,7 +1020,79 @@ impl Duration {
                   without modifying the original"]
     #[inline]
     pub fn mul_f64(self, rhs: f64) -> Duration {
-        Duration::from_secs_f64(rhs * self.as_secs_f64())
+        const DURATION_BITS: u32 = 64 + Duration::SECOND.as_nanos().ilog2() + 1;
+
+        const NUM_BITS_MANT: u32 = f64::MANTISSA_DIGITS - 1;
+        const NUM_BITS_EXP: u32 = 64 - f64::MANTISSA_DIGITS;
+
+        const SIGN_BIT: u64 = 1 << 63;
+        const MANT_IMPLIED_BIT: u64 = 1 << NUM_BITS_MANT;
+        const MANT_MASK: u64 = MANT_IMPLIED_BIT - 1;
+        const EXP_MASK: u64 = (1 << NUM_BITS_EXP) - 1;
+        const EXP_BIAS: i32 = 1 - (1 << NUM_BITS_EXP) / 2;
+
+        let rhs_raw = rhs.to_bits();
+        let mant_raw = rhs_raw & MANT_MASK;
+        let exp_raw = (rhs_raw >> NUM_BITS_MANT) & EXP_MASK;
+
+        if exp_raw == EXP_MASK {
+            panic!("result is not finite");
+        }
+
+        if self.is_zero() || rhs_raw & !SIGN_BIT == 0 {
+            // At least one factor is zero, so their product must be zero.
+            return Duration::ZERO;
+        }
+
+        let (mant, biased_exp) = if exp_raw == 0 {
+            // `rhs` is subnormal.  Renormalize by shifting so that the highest set bit is
+            // in the correct leading bit position.
+            let shift = mant_raw.leading_zeros() - NUM_BITS_EXP;
+            (mant_raw << shift, 1 - shift as i32)
+        } else {
+            // `rhs` is normal.  Set the implied leading bit.
+            (mant_raw | MANT_IMPLIED_BIT, exp_raw as i32)
+        };
+
+        // From this point on, consider `mant` to be an integer of `f64::MANTISSA_DIGITS` bits,
+        // rather than a normalized fraction; therefore adjust the exponent by `NUM_BITS_MANT`.
+        let exp = biased_exp + EXP_BIAS - NUM_BITS_MANT as i32;
+        let mant = mant as u128;
+
+        let lhs_nanos = self.as_nanos();
+
+        let nanos = if let Ok(exp) = u32::try_from(exp) {
+            // `exp` is non-negative, therefore `rhs` has no fractional part.  Mere integer
+            // multiplication is sufficient.
+            mant.checked_shl(exp).and_then(|rhs_int| lhs_nanos.checked_mul(rhs_int))
+        } else {
+            let nexp = -exp as u32;
+            if nexp > DURATION_BITS + f64::MANTISSA_DIGITS {
+                // `rhs` is such a small factor that it scales down any `Duration` to less than 0.5ns.
+                return Duration::ZERO;
+            }
+
+            let (nexp_, start) = if nexp < 128 { (nexp, 0) } else { (nexp - 128, 1) };
+
+            // Shift the mantissa across three values, in preparation for long multiplication.
+            let mut split = [0; 3];
+            split[start..start + 2].copy_from_slice(&[mant >> nexp_, mant << (128 - nexp_)]);
+
+            // Perform the long multiplication.
+            lhs_nanos.checked_mul(split[0]).and_then(|product| {
+                let (_, carry) = lhs_nanos.carrying_mul(split[2], 0);
+                let (frac, carry) = lhs_nanos.carrying_mul(split[1], carry);
+                let round_up = frac >> 127;
+                product.checked_add(carry + round_up)
+            })
+        }
+        .expect("result overflows `Duration`");
+
+        if nanos != 0 && rhs.is_sign_negative() {
+            panic!("result is negative")
+        }
+
+        Duration::from_nanos_u128(nanos)
     }
 
     /// Multiplies `Duration` by `f32`.
@@ -1033,15 +1105,15 @@ impl Duration {
     /// use std::time::Duration;
     ///
     /// let dur = Duration::new(2, 700_000_000);
-    /// assert_eq!(dur.mul_f32(3.14), Duration::new(8, 478_000_641));
+    /// assert_eq!(dur.mul_f32(3.14), Duration::new(8, 478_000_283));
     /// assert_eq!(dur.mul_f32(3.14e5), Duration::new(847_800, 0));
     /// ```
     #[stable(feature = "duration_float", since = "1.38.0")]
     #[must_use = "this returns the result of the operation, \
                   without modifying the original"]
     #[inline]
     pub fn mul_f32(self, rhs: f32) -> Duration {
-        Duration::from_secs_f32(rhs * self.as_secs_f32())
+        self.mul_f64(rhs.into())
     }
 
     /// Divides `Duration` by `f64`.
@@ -1062,7 +1134,98 @@ impl Duration {
                   without modifying the original"]
     #[inline]
     pub fn div_f64(self, rhs: f64) -> Duration {
-        Duration::from_secs_f64(self.as_secs_f64() / rhs)
+        const DURATION_BITS: u32 = 64 + Duration::SECOND.as_nanos().ilog2() + 1;
+
+        const NUM_BITS_MANT: u32 = f64::MANTISSA_DIGITS - 1;
+        const NUM_BITS_EXP: u32 = 64 - f64::MANTISSA_DIGITS;
+
+        const SIGN_BIT: u64 = 1 << 63;
+        const MANT_IMPLIED_BIT: u64 = 1 << NUM_BITS_MANT;
+        const MANT_MASK: u64 = MANT_IMPLIED_BIT - 1;
+        const EXP_MASK: u64 = (1 << NUM_BITS_EXP) - 1;
+        const EXP_BIAS: i32 = 1 - (1 << NUM_BITS_EXP) / 2;
+
+        let rhs_raw = rhs.to_bits();
+
+        if rhs.is_nan() || rhs_raw & !SIGN_BIT == 0 {
+            // `rhs` is either NaN or zero.
+            panic!("result is not finite");
+        }
+
+        let mant_raw = rhs_raw & MANT_MASK;
+        let exp_raw = (rhs_raw >> NUM_BITS_MANT) & EXP_MASK;
+
+        if self.is_zero() || exp_raw == EXP_MASK {
+            // Dividend is zero and/or divisor is infinite.
+            return Duration::ZERO;
+        }
+
+        let (mant, biased_exp) = if exp_raw == 0 {
+            // `rhs` is subnormal.  Renormalize by shifting so that the highest set bit is
+            // in the correct leading bit position.
+            let shift = mant_raw.leading_zeros() - NUM_BITS_EXP;
+            (mant_raw << shift, 1 - shift as i32)
+        } else {
+            // `rhs` is normal.  Set the implied leading bit.
+            (mant_raw | MANT_IMPLIED_BIT, exp_raw as i32)
+        };
+
+        // From this point on, consider `mant` to be an integer of `f64::MANTISSA_DIGITS` bits,
+        // rather than a normalized fraction; therefore adjust the exponent by `NUM_BITS_MANT`.
+        let exp = biased_exp + EXP_BIAS - NUM_BITS_MANT as i32;
+        let mant = mant as u128;
+
+        let lhs_nanos = self.as_nanos();
+
+        let nanos = if let Ok(exp) = u32::try_from(exp) {
+            // `exp` is non-negative, therefore `rhs` has no fractional part.  Mere integer
+            // division is sufficient.
+            (lhs_nanos / mant).unbounded_shr(exp)
+        } else {
+            // `exp` is negative, therefore `rhs` has both integer and fractional parts.
+
+            // Calculate the first 127 bits of the reciprocal, and any remainder.
+            const MAX: u128 = 1 << 127;
+            let recip_q = MAX / mant;
+            let recip_r = MAX % mant;
+
+            // Calculate the next 127 bits of the reciprocal, and shift them into the most
+            // significant bits of a u128.  Whilst the least significant bit may be incorrect, it is
+            // immaterial to the following calculation (it never funnels into `split`).
+            let recip_f = ((recip_q * recip_r) + (recip_r * recip_r) / mant) << 1;
+
+            // Offset exponent by 1 because reciprocal is calculated as 2^127/mant,
+            // rather than 2^128/mant (which obviously requires greater width than u128).
+            let nexp = (1 - exp) as u32;
+            if nexp > DURATION_BITS + f64::MANTISSA_DIGITS {
+                // `rhs` is such a small divisor that it scales up any `Duration` to overflow.
+                panic!("result overflows `Duration`");
+            }
+
+            let (nexp_, start) = if nexp < 128 { (nexp, 1) } else { (nexp - 128, 0) };
+
+            // Shift the reciprocal across three values, in preparation for long multiplication.
+            let mut split = [0; 3];
+            split[start..start + 2]
+                .copy_from_slice(&[recip_q >> (128 - nexp_), recip_q.funnel_shl(recip_f, nexp_)]);
+
+            // Perform the long multiplication.
+            lhs_nanos
+                .checked_mul(split[1])
+                .and_then(|product| {
+                    let (_, carry) = lhs_nanos.carrying_mul(split[0], 0);
+                    let (frac, carry) = lhs_nanos.carrying_mul(split[2], carry);
+                    let round_up = frac >> 127;
+                    product.checked_add(carry + round_up)
+                })
+                .expect("result overflows `Duration`")
+        };
+
+        if nanos != 0 && rhs.is_sign_negative() {
+            panic!("result is negative")
+        }
+
+        Duration::from_nanos_u128(nanos)
     }
 
     /// Divides `Duration` by `f32`.
@@ -1077,15 +1240,15 @@ impl Duration {
     /// let dur = Duration::new(2, 700_000_000);
     /// // note that due to rounding errors result is slightly
     /// // different from 0.859_872_611
-    /// assert_eq!(dur.div_f32(3.14), Duration::new(0, 859_872_580));
+    /// assert_eq!(dur.div_f32(3.14), Duration::new(0, 859_872_583));
     /// assert_eq!(dur.div_f32(3.14e5), Duration::new(0, 8_599));
     /// ```
     #[stable(feature = "duration_float", since = "1.38.0")]
     #[must_use = "this returns the result of the operation, \
                   without modifying the original"]
     #[inline]
     pub fn div_f32(self, rhs: f32) -> Duration {
-        Duration::from_secs_f32(self.as_secs_f32() / rhs)
+        self.div_f64(rhs.into())
     }
 
     /// Divides `Duration` by `Duration` and returns `f64`.

library/coretests/tests/time.rs

@@ -598,3 +598,70 @@ fn from_neg_zero() {
     assert_eq!(Duration::from_secs_f32(-0.0), Duration::ZERO);
     assert_eq!(Duration::from_secs_f64(-0.0), Duration::ZERO);
 }
+
+#[test]
+fn precise_duration_fp_mul() {
+    let d1 = Duration::from_nanos_u128(1 << 90);
+    let d2 = Duration::from_nanos_u128(2 << 90);
+    let d3 = Duration::from_nanos_u128(3 << 90);
+
+    assert_eq!(d1.mul_f32(1.0), d1);
+    assert_eq!(d1.mul_f32(2.0), d2);
+    assert_eq!(d1.mul_f32(3.0), d3);
+    assert_eq!(d2.mul_f32(1.5), d3);
+
+    let _ = Duration::MAX.mul_f32(1.0);
+}
+
+#[test]
+fn precise_duration_fp_div() {
+    let d1 = Duration::from_nanos_u128(1 << 90);
+    let d2 = Duration::from_nanos_u128(2 << 90);
+    let d3 = Duration::from_nanos_u128(3 << 90);
+
+    assert_eq!(d1.div_f32(1.0), d1);
+    assert_eq!(d2.div_f32(2.0), d1);
+    assert_eq!(d3.div_f32(3.0), d1);
+    assert_eq!(d3.div_f32(1.5), d2);
+
+    let _ = Duration::MAX.div_f32(1.0);
+}
+
+const TOO_LARGE: f64 = Duration::MAX.as_nanos() as f64;
+
+// Smallest factor is half the smallest divisor due to rounding in the opposite direction.
+const SMALLEST_DIVISOR: f64 = TOO_LARGE.recip();
+const SMALLEST_FACTOR: f64 = SMALLEST_DIVISOR / 2.0;
+
+#[test]
+fn precise_duration_fp_boundaries() {
+    const DURATION_BITS: u32 = 64 + Duration::SECOND.as_nanos().ilog2() + 1;
+    const PRECISION: u32 = DURATION_BITS - f64::MANTISSA_DIGITS;
+
+    assert_eq!(Duration::MAX.mul_f64(SMALLEST_FACTOR), Duration::NANOSECOND);
+    assert_eq!(Duration::MAX.mul_f64(SMALLEST_FACTOR.next_down()), Duration::ZERO);
+    assert_eq!(Duration::MAX.div_f64(TOO_LARGE), Duration::ZERO);
+    assert_eq!(Duration::MAX.div_f64(TOO_LARGE.next_down()), Duration::NANOSECOND);
+
+    // the following assertions pair with the two subsequent (`should_panic`) tests
+    assert!(
+        Duration::MAX - Duration::NANOSECOND.mul_f64(TOO_LARGE.next_down())
+            < Duration::from_nanos(1 << PRECISION)
+    );
+    assert!(
+        Duration::MAX - Duration::NANOSECOND.div_f64(SMALLEST_DIVISOR)
+            < Duration::from_nanos(1 << PRECISION)
+    );
+}
+
+#[test]
+#[should_panic(expected = "overflow")]
+fn precise_duration_fp_mul_overflow() {
+    let _ = Duration::NANOSECOND.mul_f64(TOO_LARGE);
+}
+
+#[test]
+#[should_panic(expected = "overflow")]
+fn precise_duration_fp_div_overflow() {
+    let _ = Duration::NANOSECOND.div_f64(SMALLEST_DIVISOR.next_down());
+}