Improve accuracy of asinh and acosh

Author: mwillsey

library/std/src/f32.rs

@@ -880,7 +880,9 @@ impl f32 {
     #[stable(feature = "rust1", since = "1.0.0")]
     #[inline]
     pub fn asinh(self) -> f32 {
-        (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self)
+        let ax = self.abs();
+        let ix = 1.0 / ax;
+        (ax + (ax / (Self::hypot(1.0, ix) + ix))).ln_1p().copysign(self)
     }
 
     /// Inverse hyperbolic cosine function.
@@ -900,7 +902,11 @@ impl f32 {
     #[stable(feature = "rust1", since = "1.0.0")]
     #[inline]
     pub fn acosh(self) -> f32 {
-        if self < 1.0 { Self::NAN } else { (self + ((self * self) - 1.0).sqrt()).ln() }
+        if self < 1.0 {
+            Self::NAN
+        } else {
+            (self + ((self - 1.0).sqrt() * (self + 1.0).sqrt())).ln()
+        }
     }
 
     /// Inverse hyperbolic tangent function.

library/std/src/f32/tests.rs

@@ -587,6 +587,11 @@ fn test_asinh() {
     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 for low accuracy from issue 104548
+    assert_approx_eq!(60.0f32, 60.0f32.sinh().asinh());
+    // mul needed for approximate comparison to be meaningful
+    assert_approx_eq!(1.0f32, 1e-15f32.sinh().asinh() * 1e15f32);
 }
 
 #[test]
@@ -602,6 +607,9 @@ fn test_acosh() {
     assert!(nan.acosh().is_nan());
     assert_approx_eq!(2.0f32.acosh(), 1.31695789692481670862504634730796844f32);
     assert_approx_eq!(3.0f32.acosh(), 1.76274717403908605046521864995958461f32);
+
+    // test for low accuracy from issue 104548
+    assert_approx_eq!(60.0f32, 60.0f32.cosh().acosh());
 }
 
 #[test]

library/std/src/f64.rs

@@ -882,7 +882,9 @@ impl f64 {
     #[stable(feature = "rust1", since = "1.0.0")]
     #[inline]
     pub fn asinh(self) -> f64 {
-        (self.abs() + ((self * self) + 1.0).sqrt()).ln().copysign(self)
+        let ax = self.abs();
+        let ix = 1.0 / ax;
+        (ax + (ax / (Self::hypot(1.0, ix) + ix))).ln_1p().copysign(self)
     }
 
     /// Inverse hyperbolic cosine function.
@@ -902,7 +904,11 @@ impl f64 {
     #[stable(feature = "rust1", since = "1.0.0")]
     #[inline]
     pub fn acosh(self) -> f64 {
-        if self < 1.0 { Self::NAN } else { (self + ((self * self) - 1.0).sqrt()).ln() }
+        if self < 1.0 {
+            Self::NAN
+        } else {
+            (self + ((self - 1.0).sqrt() * (self + 1.0).sqrt())).ln()
+        }
     }
 
     /// Inverse hyperbolic tangent function.

library/std/src/f64/tests.rs

@@ -575,6 +575,11 @@ fn test_asinh() {
     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 for low accuracy from issue 104548
+    assert_approx_eq!(60.0f64, 60.0f64.sinh().asinh());
+    // mul needed for approximate comparison to be meaningful
+    assert_approx_eq!(1.0f64, 1e-15f64.sinh().asinh() * 1e15f64);
 }
 
 #[test]
@@ -590,6 +595,9 @@ fn test_acosh() {
     assert!(nan.acosh().is_nan());
     assert_approx_eq!(2.0f64.acosh(), 1.31695789692481670862504634730796844f64);
     assert_approx_eq!(3.0f64.acosh(), 1.76274717403908605046521864995958461f64);
+
+    // test for low accuracy from issue 104548
+    assert_approx_eq!(60.0f64, 60.0f64.cosh().acosh());
 }
 
 #[test]