diff options
| author | Szabolcs Nagy <nsz@port70.net> | 2023-08-18 23:16:00 +0200 |
|---|---|---|
| committer | Rich Felker <dalias@aerifal.cx> | 2023-08-19 14:07:44 -0400 |
| commit | 6d10102709df4bc966d2846c1c45cd667e5048e5 (patch) | |
| tree | 6fd927698010a49a910bfed6923f882b909015ec /src | |
| parent | 83b858f83b658bd34eca5d8ad4d145f673ae7e5e (diff) | |
| download | musl-6d10102709df4bc966d2846c1c45cd667e5048e5.tar.gz | |
math: fix ld80 acoshl(x) for x < 0
acosh(x) is nan for x < 1, but x < 0 cases were not handled specially
and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.:
acoshl(-0x1p20) returned -inf,
acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L,
fixed by checking the sign bit and handling it specially.
reported by Paul Zimmermann.
Diffstat (limited to 'src')
| -rw-r--r-- | src/math/acoshl.c | 10 |
1 files changed, 7 insertions, 3 deletions
diff --git a/src/math/acoshl.c b/src/math/acoshl.c index 8d4b43f6..943cec17 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -10,14 +10,18 @@ long double acoshl(long double x) long double acoshl(long double x) { union ldshape u = {x}; - int e = u.i.se & 0x7fff; + int e = u.i.se; if (e < 0x3fff + 1) - /* |x| < 2, invalid if x < 1 or nan */ + /* 0 <= x < 2, invalid if x < 1 */ return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); if (e < 0x3fff + 32) - /* |x| < 0x1p32 */ + /* 2 <= x < 0x1p32 */ return logl(2*x - 1/(x+sqrtl(x*x-1))); + if (e & 0x8000) + /* x < 0 or x = -0, invalid */ + return (x - x) / (x - x); + /* 0x1p32 <= x or nan */ return logl(x) + 0.693147180559945309417232121458176568L; } #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |