Fix asinh of negative values

[Ralith]

Rust's current implementation of asinh has large errors in its negative range. These are (mostly) not numerical, but rather seem due to an incorrect implementation. This appears to be due to avoidable catastrophic cancellation.
Playground before/after.
glibc uses abs here.

Many thanks to @danieldeankon for finding this weird behavior, @jebrosen for diagnosing it, and @toasteater for identifying the probable implementation error!

[rust-highfive]

r? @shepmaster

(rust_highfive has picked a reviewer for you, use r? to override)

[ghost]

To think about it after the fact, this probably is a numerical issue, though. self + ((self * self) + 1.0).sqrt() is supposed to produce a very small value whose .ln() is equal to -(self.abs() + ((self * self) + 1.0).sqrt()).ln(). In practice, however, once the absolute value of self goes above 2^26 (67108864, compare to 67452095 in the original example), self * self will be around 2^53, and the +1.0 will no longer work. That makes self + ((self * self) + 1.0).sqrt() dangerously close (or equal to) zero, leading to the precision problems and -infs in the original example.

The .abs() addition is a transformation that works because asinh is odd, and additions are much more permissive when there are no negative numbers involved.

[shepmaster]

I don't have the numerical chops to evaluate this, so I'll be reassigning, but it seems like a really bad idea to "fix" this without (a) even stating what the problem is in the commit message (links to external sites are great, but ultimately very fragile) (b) adding any kind of automated test to prevent regressions.

r? @cramertj

[Ralith]

Good points. Added discussion of the numerical issue and rationale for the correctness of the solution to the commit, and threw in a basic regression test.

[Ralith]

Ping?

[Mark-Simulacrum]

Hm, I'm not sure who a good reviewer here would be.

Cc @rust-lang/libs - do we have anyone?

[dtolnay]

I can review this kind of thing.

Asinh is odd, and it's clear we're getting better results in the positive range, so I'm happy to take this improvement. libclc does the abs there too.

// Arguments such that 4.0 <= abs(x) <= 1/sqrt(epsilon) are
// approximated by asinhf(x) = ln(abs(x) + sqrt(x*x+1))
// with the sign of x (see Abramowitz and Stegun 4.6.20)

Following up on this, it would be good to look at an appropriately licensed high quality implementation and consider whether we should pick up what they do on some of the other intervals, namely close to zero and very far from zero where ln(abs(x) + sqrt(x*x + 1)) still isn't good.

[dtolnay]

@bors r+

[bors]

📌 Commit 730f736 has been approved by dtolnay

[dtolnay]

@bors r+

[bors]

📌 Commit 35a2915 has been approved by dtolnay

[RalfJung]

@bors rollup