From 343e091465e81589a164adba621f331201746dc6 Mon Sep 17 00:00:00 2001 From: Travis Cross Date: Mon, 2 Dec 2024 13:22:48 +0000 Subject: [PATCH] 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. --- library/std/src/f128.rs | 67 +++++++++++++++++++++++++++++++++++++---- library/std/src/f16.rs | 63 +++++++++++++++++++++++++++++++++++--- library/std/src/f32.rs | 63 ++++++++++++++++++++++++++++++++++---- library/std/src/f64.rs | 63 ++++++++++++++++++++++++++++++++++---- 4 files changed, 234 insertions(+), 22 deletions(-) diff --git a/library/std/src/f128.rs b/library/std/src/f128.rs index e93e915159e40..56d340ab8d5fb 100644 --- a/library/std/src/f128.rs +++ b/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 { + 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)`. diff --git a/library/std/src/f16.rs b/library/std/src/f16.rs index 5b7fcaa28e064..b56b8d8359389 100644 --- a/library/std/src/f16.rs +++ b/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 { + 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)`. diff --git a/library/std/src/f32.rs b/library/std/src/f32.rs index 7cb285bbff5f7..ab63bec01eb69 100644 --- a/library/std/src/f32.rs +++ b/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 { + 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)`. diff --git a/library/std/src/f64.rs b/library/std/src/f64.rs index 47163c272de32..92564df35ea70 100644 --- a/library/std/src/f64.rs +++ b/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 { + 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)`.