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

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)`.