Skip to content

Arc-hyperbolic-trigonometric functions are very inaccurate #104548

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
pavpanchekha opened this issue Nov 17, 2022 · 4 comments · Fixed by #104553
Closed

Arc-hyperbolic-trigonometric functions are very inaccurate #104548

pavpanchekha opened this issue Nov 17, 2022 · 4 comments · Fixed by #104553
Labels
C-bug Category: This is a bug.

Comments

@pavpanchekha
Copy link
Contributor

pavpanchekha commented Nov 17, 2022

The functions asinh and acosh for both f32 and f64 are quite inaccurate, both for very large inputs (because x * x overflows) and, for asinh, for very small inputs (due to kinda-complex cancellation).

I tried this code:

fn main() {
    println!("{}", (600.0 as f64).sinh().asinh());
    println!("{}", (1e-16 as f64).sinh().asinh());
    println!("{}", (600.0 as f64).cosh().acosh());
    
    println!("{}", (60.0 as f32).sinh().asinh());
    println!("{}", (1e-8 as f32).sinh().asinh());
    println!("{}", (60.0 as f32).cosh().acosh());
}

This code evaluates asinh(sinh(x)) and acosh(cosh(x)). Naturally I expected the output to be x, but instead the output is inf and 0, indicating overflow and cancellation.

Meta

I ran this in the Rust Playground on both Stable (1.65.0) and Nightly (2022-11-16 e9493d63c2a57b91556d).

@pavpanchekha pavpanchekha added the C-bug Category: This is a bug. label Nov 17, 2022
@pavpanchekha
Copy link
Contributor Author

pavpanchekha commented Nov 17, 2022

The issue was originally discovered by @finnbear in FPBench/FPBench#118.

I'm testing solutions that don't rely on cmath in herbie-fp/herbie#514.

@pavpanchekha
Copy link
Contributor Author

In the two PRs above the following implementations are proposed:

fn asinh(x: f64) -> f64 {
	let ax = x.abs();
	let ix = 1.0 / ax;
	((ax + (ax / (f64::hypot(1.0, ix) + ix)))).ln_1p().copysign(x)
}

fn acosh(x: f64) -> f64 {
	if x < 1.0 { f64::NAN } else { (x + ((x - 1.0).sqrt() * (x + 1.0).sqrt())).ln() }
}

I'll file a PR once I get time to figure out how to download + build + test Rust, probably this or next week.

@pavpanchekha
Copy link
Contributor Author

Here is the accuracy plot from Herbie for f64::asinh, red being the current Rust expression and green being the proposed fix:

image

And here's f64::acosh:

image

The f32 versions look similar but the problems occur at different (smaller) inputs.

@finnbear
Copy link
Contributor

finnbear commented Nov 17, 2022

I'll file a PR once I get time to figure out how to download + build + test Rust, probably this or next week.

I already have a Rust fork and the required tools to test it so, if it's ok with you, I could write the PR using your proposed implementation.

Looks like #104553 beat me to it 😃

@bors bors closed this as completed in 45bfb1c Nov 19, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
C-bug Category: This is a bug.
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants