math: fix acosh for x < 1 (including negative inputs) #35616
arshidkv12/zig:acosh-1 into master
I don't think this is the fix we want here.
Specifically, I expect the x < 1.0 case to be handled in the first branch:
return math.log1p(x - 1 + @sqrt((x - 1) * (x - 1) + 2 * (x - 1)));
This handles the cases implicitly by the behaviour of log1p and sqrt. Either the argument to sqrt will be negative (-1 < x < 1), or the argument to log1p will be negative (x < -1) ), as x will grow faster relative to sqrt argument in all other cases. Both of which will propagate a nan through as the acosh return value.
If you run these tests against master, we can see that f32 already produces the correct values according to the tests, suggesting that only acosh64 is an issue.
Based on your debugging, perhaps look at why the second branch is being taken for these negative instead of the first (which should handle the domain issue). Notably the exponent check differs a little.
See also upstream musl which this code was derived from: https://git.musl-libc.org/cgit/musl/tree/src/math/acosh.c
@tiehuis wrote in #35616 (comment):
I don't think this is the fix we want here.
Specifically, I expect the
x < 1.0case to be handled in the first branch:return math.log1p(x - 1 + @sqrt((x - 1) * (x - 1) + 2 * (x - 1)));This handles the cases implicitly by the behaviour of log1p and sqrt. Either the argument to sqrt will be negative (
-1 < x < 1), or the argument tolog1pwill be negative (x < -1) ), as x will grow faster relative tosqrtargument in all other cases. Both of which will propagate a nan through as the acosh return value.If you run these tests against master, we can see that
f32already produces the correct values according to the tests, suggesting that onlyacosh64is an issue.Based on your debugging, perhaps look at why the second branch is being taken for these negative instead of the first (which should handle the domain issue). Notably the exponent check differs a little.
See also upstream musl which this code was derived from: https://git.musl-libc.org/cgit/musl/tree/src/math/acosh.c
Yes, the acosh32 is working properly. But I think, musl c code also not working.
#include <stdio.h>#include <math.h>#include <stdint.h>
double acosh(double x)
{
union {double f; uint64_t i;} u = {.f = x};
unsigned e = u.i >> 52 & 0x7ff;
/* x < 1 domain error is handled in the called functions */
if (e < 0x3ff + 1)
/* |x| < 2, up to 2ulp error in [1,1.125] */
return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1)));
if (e < 0x3ff + 26)
/* |x| < 0x1p26 */
return log(2*x - 1/(x+sqrt(x*x-1)));
/* |x| >= 0x1p26 or nan */
return log(x) + 0.693147180559945309417232121458176568;
}
int main(){
double x = -1.9e4;
double z = acosh(x);
printf("out x %f \n", z);
double y = -2e4;
double k = acosh(x);
printf("out x %f \n", k);
}
Out
gcc z.c && ./a.out
out x -7.288321
out x -7.288321
Note: x < 1 domain error is handled in the called functions.
Given we are have some inherent issues with the implementation's accuracy itself, I would be keen to replace the implementation and check that against the libc-test you've got set up. We would feel a bit more confident in the results instead of just handling the 0 case here. This would also help your other pr and we could enable the tests in that PR.
The FreeBSD implementation suggested in the issue is worth a first look.
Thank you. I changed it based on FreeBSD.
View command line instructions
Checkout
From your project repository, check out a new branch and test the changes.No due date set.
No dependencies set.
Deleting a branch is permanent. Although the deleted branch may continue to exist for a short time before it actually gets removed, it CANNOT be undone in most cases. Continue?