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.

1 files changed, 3 insertions, 1 deletions

diff --git a/src/libstd/f64.rs b/src/libstd/f64.rs
index fe64d27b1ef..c033198f021 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]