From 730f7366bba5abdf5ae0c2f1222795e40d48f90c Mon Sep 17 00:00:00 2001 From: Benjamin Saunders Date: Fri, 22 May 2020 20:49:20 -0700 Subject: [PATCH 1/2] Fix asinh of negative values When `x` has large magnitude, `x + ((x * x) + 1.0).sqrt()` approaches `x + x.abs()`. For negative values of `x`, this leads to catastrophic cancellation, resulting in large errors or even 0 being passed to `ln`, producing incorrect results including `-inf`. Becuase asinh is an odd function, i.e. -asinh(x) = asinh(-x) for all x, we can avoid the catastrophic cancellation and obtain correct results by taking the absolute value of `self` for the first term. `self * self` is always positive, so in effect this gives us `x.abs().asinh().copysign(x)` which as discussed above is algebraically equivalent, but is much more accurate. --- src/libstd/f32.rs | 4 +++- src/libstd/f64.rs | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/libstd/f32.rs b/src/libstd/f32.rs index 8e743ace99bfb..ab468f42b3bd0 100644 --- a/src/libstd/f32.rs +++ b/src/libstd/f32.rs @@ -835,7 +835,7 @@ impl f32 { if self == Self::NEG_INFINITY { Self::NEG_INFINITY } else { - (self + ((self * self) + 1.0).sqrt()).ln().copysign(self) + (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) } } @@ -1414,6 +1414,8 @@ mod tests { assert!((-0.0f32).asinh().is_sign_negative()); // issue 63271 assert_approx_eq!(2.0f32.asinh(), 1.443635475178810342493276740273105f32); assert_approx_eq!((-2.0f32).asinh(), -1.443635475178810342493276740273105f32); + // regression test for the catastrophic cancellation fixed in 72486 + assert_approx_eq!((-3000.0f32).asinh(), -8.699514775987968673236893537700647f32); } #[test] diff --git a/src/libstd/f64.rs b/src/libstd/f64.rs index fe64d27b1efc8..c033198f021c2 100644 --- a/src/libstd/f64.rs +++ b/src/libstd/f64.rs @@ -837,7 +837,7 @@ impl f64 { if self == Self::NEG_INFINITY { Self::NEG_INFINITY } else { - (self + ((self * self) + 1.0).sqrt()).ln().copysign(self) + (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) } } @@ -1443,6 +1443,8 @@ mod tests { // issue 63271 assert_approx_eq!(2.0f64.asinh(), 1.443635475178810342493276740273105f64); assert_approx_eq!((-2.0f64).asinh(), -1.443635475178810342493276740273105f64); + // regression test for the catastrophic cancellation fixed in 72486 + assert_approx_eq!((-67452098.07139316f64).asinh(), -18.72007542627454439398548429400083); } #[test] From 35a2915bf3f24edbbee2d6a9c8f2318304b99df3 Mon Sep 17 00:00:00 2001 From: Benjamin Saunders Date: Thu, 18 Jun 2020 18:52:26 -0700 Subject: [PATCH 2/2] Remove now-redundant branch --- src/libstd/f32.rs | 6 +----- src/libstd/f64.rs | 6 +----- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/src/libstd/f32.rs b/src/libstd/f32.rs index ab468f42b3bd0..7899c422da62d 100644 --- a/src/libstd/f32.rs +++ b/src/libstd/f32.rs @@ -832,11 +832,7 @@ impl f32 { #[stable(feature = "rust1", since = "1.0.0")] #[inline] pub fn asinh(self) -> f32 { - if self == Self::NEG_INFINITY { - Self::NEG_INFINITY - } else { - (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) - } + (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) } /// Inverse hyperbolic cosine function. diff --git a/src/libstd/f64.rs b/src/libstd/f64.rs index c033198f021c2..c663c5ffc9beb 100644 --- a/src/libstd/f64.rs +++ b/src/libstd/f64.rs @@ -834,11 +834,7 @@ impl f64 { #[stable(feature = "rust1", since = "1.0.0")] #[inline] pub fn asinh(self) -> f64 { - if self == Self::NEG_INFINITY { - Self::NEG_INFINITY - } else { - (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) - } + (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self) } /// Inverse hyperbolic cosine function.