diff options
| author | Benjamin Saunders <ben.e.saunders@gmail.com> | 2020-05-22 20:49:20 -0700 |
|---|---|---|
| committer | Benjamin Saunders <ben.e.saunders@gmail.com> | 2020-05-24 00:03:35 -0700 |
| commit | 730f7366bba5abdf5ae0c2f1222795e40d48f90c (patch) | |
| tree | 84536380e5fc3796d268a5f3e7ed9c0e770bd9de | |
| parent | 215f2d3294b08dbdcf8f7d40de21ef1e7eae0a2d (diff) | |
| download | rust-730f7366bba5abdf5ae0c2f1222795e40d48f90c.tar.gz rust-730f7366bba5abdf5ae0c2f1222795e40d48f90c.zip | |
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.
| -rw-r--r-- | src/libstd/f32.rs | 4 | ||||
| -rw-r--r-- | 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 8e743ace99b..ab468f42b3b 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 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] |
