Draft: Fix {f16,f32,f64,f128}::div_euclid

The current implementation of `div_euclid` for floating point numbers
violates the invariant, stated in the documentation, that:

```rust
a.rem_euclid(b) ~= a.div_euclid(b).mul_add(-b, a)
```

This fixes that problem.

When the magnitude of the exact quotient is greater than or equal to
the maximum integer that can be represented precisely in the type,
that invariant is not generally possible to uphold -- and without
control over the rounding mode it would be difficult to calculate in
any case -- so we return `NaN` in those instances.

Author: traviscross

library/std/src/f128.rs

@@ -236,15 +236,19 @@ impl f128 {
     /// Calculates Euclidean division, the matching method for `rem_euclid`.
     ///
     /// This computes the integer `n` such that
-    /// `self = n * rhs + self.rem_euclid(rhs)`.
-    /// In other words, the result is `self / rhs` rounded to the integer `n`
-    /// such that `self >= n * rhs`.
+    /// `self.rem_euclid(rhs) = self.div_euclid(rhs).mul_add(-rhs, self)`.
+    /// In other words, the result is `self / rhs` rounded to the largest
+    /// integer `n` such that `self >= n * rhs`.
     ///
     /// # Precision
     ///
     /// The result of this operation is guaranteed to be the rounded
     /// infinite-precision result.
     ///
+    /// If the magnitude of the exact quotient `self / rhs` is greater than or
+    /// equal to the maximum integer that can be represented precisely with this
+    /// type, then `NaN` may be returned.
+    ///
     /// # Examples
     ///
     /// ```
@@ -259,16 +263,67 @@ impl f128 {
     /// assert_eq!((-a).div_euclid(-b), 2.0); // -7.0 >= -4.0 * 2.0
     /// # }
     /// ```
+    ///
+    /// ```
+    /// #![feature(f128)]
+    /// # #[cfg(reliable_f128_math)] {
+    ///
+    /// let a: f32 = 11.0;
+    /// let b = 1.1;
+    /// assert_eq!(a.div_euclid(b), 9.0);
+    /// assert_eq!(a.rem_euclid(b), a.div_euclid(b).mul_add(-b, a));
+    /// # }
+    /// ```
     #[inline]
     #[rustc_allow_incoherent_impl]
     #[unstable(feature = "f128", issue = "116909")]
     #[must_use = "method returns a new number and does not mutate the original value"]
     pub fn div_euclid(self, rhs: f128) -> f128 {
-        let q = (self / rhs).trunc();
+        use core::ops::ControlFlow;
+
+        fn xor_sign(x: f128, y: f128, v: f128) -> f128 {
+            match () {
+                _ if x.is_sign_positive() == y.is_sign_positive() => v,
+                _ => -v,
+            }
+        }
+
+        fn div_special(x: f128, y: f128) -> ControlFlow<f128> {
+            ControlFlow::Break(match (x, y) {
+                (x, y) if x.is_nan() || y.is_nan() => f128::NAN,
+                (x, y) if x == 0.0 && y == 0.0 => f128::NAN,
+                (x, y) if x == 0.0 => xor_sign(x, y, 0.0),
+                (x, y) if y == 0.0 => xor_sign(x, y, f128::INFINITY),
+                (x, y) if x.is_infinite() && y.is_infinite() => f128::NAN,
+                (x, y) if x.is_infinite() => xor_sign(x, y, f128::INFINITY),
+                (x, y) if y.is_infinite() => xor_sign(x, y, 0.0),
+                _ => return ControlFlow::Continue(()),
+            })
+        }
+
+        fn div_trunc(x: f128, y: f128) -> f128 {
+            const MAX_EXACT_Q: f128 = ((1u128 << f128::MANTISSA_DIGITS) - 1) as f128;
+            if let ControlFlow::Break(q) = div_special(x, y) {
+                return q;
+            }
+            let sign = xor_sign(x, y, 1.0);
+            let (x, y) = (x.abs(), y.abs());
+            let q = x / y;
+            if q > MAX_EXACT_Q {
+                return f128::NAN;
+            }
+            let qt = q.trunc();
+            if qt.mul_add(-y, x).is_sign_negative() { qt - 1.0 } else { qt }.copysign(sign)
+        }
+
+        if let ControlFlow::Break(q) = div_special(self, rhs) {
+            return q;
+        }
+        let qt = div_trunc(self, rhs);
         if self % rhs < 0.0 {
-            return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
+            return if rhs > 0.0 { qt - 1.0 } else { qt + 1.0 };
         }
-        q
+        qt
     }
 
     /// Calculates the least nonnegative remainder of `self (mod rhs)`.

library/std/src/f16.rs

@@ -236,7 +236,7 @@ impl f16 {
     /// Calculates Euclidean division, the matching method for `rem_euclid`.
     ///
     /// This computes the integer `n` such that
-    /// `self = n * rhs + self.rem_euclid(rhs)`.
+    /// `self.rem_euclid(rhs) = self.div_euclid(rhs).mul_add(-rhs, self)`.
     /// In other words, the result is `self / rhs` rounded to the integer `n`
     /// such that `self >= n * rhs`.
     ///
@@ -245,6 +245,10 @@ impl f16 {
     /// The result of this operation is guaranteed to be the rounded
     /// infinite-precision result.
     ///
+    /// If the magnitude of the exact quotient `self / rhs` is greater than or
+    /// equal to the maximum integer that can be represented precisely with this
+    /// type, then `NaN` may be returned.
+    ///
     /// # Examples
     ///
     /// ```
@@ -259,16 +263,67 @@ impl f16 {
     /// assert_eq!((-a).div_euclid(-b), 2.0); // -7.0 >= -4.0 * 2.0
     /// # }
     /// ```
+    ///
+    /// ```
+    /// #![feature(f16)]
+    /// # #[cfg(reliable_f16_math)] {
+    ///
+    /// let a: f16 = 11.0;
+    /// let b = 1.1;
+    /// assert_eq!(a.div_euclid(b), 9.0);
+    /// assert_eq!(a.rem_euclid(b), a.div_euclid(b).mul_add(-b, a));
+    /// # }
+    /// ```
     #[inline]
     #[rustc_allow_incoherent_impl]
     #[unstable(feature = "f16", issue = "116909")]
     #[must_use = "method returns a new number and does not mutate the original value"]
     pub fn div_euclid(self, rhs: f16) -> f16 {
-        let q = (self / rhs).trunc();
+        use core::ops::ControlFlow;
+
+        fn xor_sign(x: f16, y: f16, v: f16) -> f16 {
+            match () {
+                _ if x.is_sign_positive() == y.is_sign_positive() => v,
+                _ => -v,
+            }
+        }
+
+        fn div_special(x: f16, y: f16) -> ControlFlow<f16> {
+            ControlFlow::Break(match (x, y) {
+                (x, y) if x.is_nan() || y.is_nan() => f16::NAN,
+                (x, y) if x == 0.0 && y == 0.0 => f16::NAN,
+                (x, y) if x == 0.0 => xor_sign(x, y, 0.0),
+                (x, y) if y == 0.0 => xor_sign(x, y, f16::INFINITY),
+                (x, y) if x.is_infinite() && y.is_infinite() => f16::NAN,
+                (x, y) if x.is_infinite() => xor_sign(x, y, f16::INFINITY),
+                (x, y) if y.is_infinite() => xor_sign(x, y, 0.0),
+                _ => return ControlFlow::Continue(()),
+            })
+        }
+
+        fn div_trunc(x: f16, y: f16) -> f16 {
+            const MAX_EXACT_Q: f16 = ((1u64 << f16::MANTISSA_DIGITS) - 1) as f16;
+            if let ControlFlow::Break(q) = div_special(x, y) {
+                return q;
+            }
+            let sign = xor_sign(x, y, 1.0);
+            let (x, y) = (x.abs(), y.abs());
+            let q = x / y;
+            if q > MAX_EXACT_Q {
+                return f16::NAN;
+            }
+            let qt = q.trunc();
+            if qt.mul_add(-y, x).is_sign_negative() { qt - 1.0 } else { qt }.copysign(sign)
+        }
+
+        if let ControlFlow::Break(q) = div_special(self, rhs) {
+            return q;
+        }
+        let qt = div_trunc(self, rhs);
         if self % rhs < 0.0 {
-            return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
+            return if rhs > 0.0 { qt - 1.0 } else { qt + 1.0 };
         }
-        q
+        qt
     }
 
     /// Calculates the least nonnegative remainder of `self (mod rhs)`.

library/std/src/f32.rs

@@ -220,15 +220,19 @@ impl f32 {
     /// Calculates Euclidean division, the matching method for `rem_euclid`.
     ///
     /// This computes the integer `n` such that
-    /// `self = n * rhs + self.rem_euclid(rhs)`.
-    /// In other words, the result is `self / rhs` rounded to the integer `n`
-    /// such that `self >= n * rhs`.
+    /// `self.rem_euclid(rhs) = self.div_euclid(rhs).mul_add(-rhs, self)`.
+    /// In other words, the result is `self / rhs` rounded to the largest
+    /// integer `n` such that `self >= n * rhs`.
     ///
     /// # Precision
     ///
     /// The result of this operation is guaranteed to be the rounded
     /// infinite-precision result.
     ///
+    /// If the magnitude of the exact quotient `self / rhs` is greater than or
+    /// equal to the maximum integer that can be represented precisely with this
+    /// type, then `NaN` may be returned.
+    ///
     /// # Examples
     ///
     /// ```
@@ -239,16 +243,63 @@ impl f32 {
     /// assert_eq!(a.div_euclid(-b), -1.0); // 7.0 >= -4.0 * -1.0
     /// assert_eq!((-a).div_euclid(-b), 2.0); // -7.0 >= -4.0 * 2.0
     /// ```
+    ///
+    /// ```
+    /// let a: f32 = 11.0;
+    /// let b = 1.1;
+    /// assert_eq!(a.div_euclid(b), 9.0);
+    /// assert_eq!(a.rem_euclid(b), a.div_euclid(b).mul_add(-b, a));
+    /// ```
     #[rustc_allow_incoherent_impl]
     #[must_use = "method returns a new number and does not mutate the original value"]
     #[inline]
     #[stable(feature = "euclidean_division", since = "1.38.0")]
     pub fn div_euclid(self, rhs: f32) -> f32 {
-        let q = (self / rhs).trunc();
+        use core::ops::ControlFlow;
+
+        fn xor_sign(x: f32, y: f32, v: f32) -> f32 {
+            match () {
+                _ if x.is_sign_positive() == y.is_sign_positive() => v,
+                _ => -v,
+            }
+        }
+
+        fn div_special(x: f32, y: f32) -> ControlFlow<f32> {
+            ControlFlow::Break(match (x, y) {
+                (x, y) if x.is_nan() || y.is_nan() => f32::NAN,
+                (x, y) if x == 0.0 && y == 0.0 => f32::NAN,
+                (x, y) if x == 0.0 => xor_sign(x, y, 0.0),
+                (x, y) if y == 0.0 => xor_sign(x, y, f32::INFINITY),
+                (x, y) if x.is_infinite() && y.is_infinite() => f32::NAN,
+                (x, y) if x.is_infinite() => xor_sign(x, y, f32::INFINITY),
+                (x, y) if y.is_infinite() => xor_sign(x, y, 0.0),
+                _ => return ControlFlow::Continue(()),
+            })
+        }
+
+        fn div_trunc(x: f32, y: f32) -> f32 {
+            const MAX_EXACT_Q: f32 = ((1u64 << f32::MANTISSA_DIGITS) - 1) as f32;
+            if let ControlFlow::Break(q) = div_special(x, y) {
+                return q;
+            }
+            let sign = xor_sign(x, y, 1.0);
+            let (x, y) = (x.abs(), y.abs());
+            let q = x / y;
+            if q > MAX_EXACT_Q {
+                return f32::NAN;
+            }
+            let qt = q.trunc();
+            if qt.mul_add(-y, x).is_sign_negative() { qt - 1.0 } else { qt }.copysign(sign)
+        }
+
+        if let ControlFlow::Break(q) = div_special(self, rhs) {
+            return q;
+        }
+        let qt = div_trunc(self, rhs);
         if self % rhs < 0.0 {
-            return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
+            return if rhs > 0.0 { qt - 1.0 } else { qt + 1.0 };
         }
-        q
+        qt
     }
 
     /// Calculates the least nonnegative remainder of `self (mod rhs)`.

library/std/src/f64.rs

@@ -220,15 +220,19 @@ impl f64 {
     /// Calculates Euclidean division, the matching method for `rem_euclid`.
     ///
     /// This computes the integer `n` such that
-    /// `self = n * rhs + self.rem_euclid(rhs)`.
-    /// In other words, the result is `self / rhs` rounded to the integer `n`
-    /// such that `self >= n * rhs`.
+    /// `self.rem_euclid(rhs) = self.div_euclid(rhs).mul_add(-rhs, self)`.
+    /// In other words, the result is `self / rhs` rounded to the largest
+    /// integer `n` such that `self >= n * rhs`.
     ///
     /// # Precision
     ///
     /// The result of this operation is guaranteed to be the rounded
     /// infinite-precision result.
     ///
+    /// If the magnitude of the exact quotient `self / rhs` is greater than or
+    /// equal to the maximum integer that can be represented precisely with this
+    /// type, then `NaN` may be returned.
+    ///
     /// # Examples
     ///
     /// ```
@@ -239,16 +243,63 @@ impl f64 {
     /// assert_eq!(a.div_euclid(-b), -1.0); // 7.0 >= -4.0 * -1.0
     /// assert_eq!((-a).div_euclid(-b), 2.0); // -7.0 >= -4.0 * 2.0
     /// ```
+    ///
+    /// ```
+    /// let a: f64 = 11.0;
+    /// let b = 1.1;
+    /// assert_eq!(a.div_euclid(b), 9.0);
+    /// assert_eq!(a.rem_euclid(b), a.div_euclid(b).mul_add(-b, a));
+    /// ```
     #[rustc_allow_incoherent_impl]
     #[must_use = "method returns a new number and does not mutate the original value"]
     #[inline]
     #[stable(feature = "euclidean_division", since = "1.38.0")]
     pub fn div_euclid(self, rhs: f64) -> f64 {
-        let q = (self / rhs).trunc();
+        use core::ops::ControlFlow;
+
+        fn xor_sign(x: f64, y: f64, v: f64) -> f64 {
+            match () {
+                _ if x.is_sign_positive() == y.is_sign_positive() => v,
+                _ => -v,
+            }
+        }
+
+        fn div_special(x: f64, y: f64) -> ControlFlow<f64> {
+            ControlFlow::Break(match (x, y) {
+                (x, y) if x.is_nan() || y.is_nan() => f64::NAN,
+                (x, y) if x == 0.0 && y == 0.0 => f64::NAN,
+                (x, y) if x == 0.0 => xor_sign(x, y, 0.0),
+                (x, y) if y == 0.0 => xor_sign(x, y, f64::INFINITY),
+                (x, y) if x.is_infinite() && y.is_infinite() => f64::NAN,
+                (x, y) if x.is_infinite() => xor_sign(x, y, f64::INFINITY),
+                (x, y) if y.is_infinite() => xor_sign(x, y, 0.0),
+                _ => return ControlFlow::Continue(()),
+            })
+        }
+
+        fn div_trunc(x: f64, y: f64) -> f64 {
+            const MAX_EXACT_Q: f64 = ((1u64 << f64::MANTISSA_DIGITS) - 1) as f64;
+            if let ControlFlow::Break(q) = div_special(x, y) {
+                return q;
+            }
+            let sign = xor_sign(x, y, 1.0);
+            let (x, y) = (x.abs(), y.abs());
+            let q = x / y;
+            if q > MAX_EXACT_Q {
+                return f64::NAN;
+            }
+            let qt = q.trunc();
+            if qt.mul_add(-y, x).is_sign_negative() { qt - 1.0 } else { qt }.copysign(sign)
+        }
+
+        if let ControlFlow::Break(q) = div_special(self, rhs) {
+            return q;
+        }
+        let qt = div_trunc(self, rhs);
         if self % rhs < 0.0 {
-            return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
+            return if rhs > 0.0 { qt - 1.0 } else { qt + 1.0 };
         }
-        q
+        qt
     }
 
     /// Calculates the least nonnegative remainder of `self (mod rhs)`.