[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <cf553cb6-cedc-72e1-1d19-4aecdeeedc3@esi.com.au>
Date: Fri, 31 Jul 2026 10:01:38 +1000 (AEST)
From: Damian McGuckin <damianm@....com.au>
To: musl@...ts.openwall.com
cc: Paul Zimmermann <Paul.Zimmermann@...ia.fr>
Subject: Re: issue in acosh
On Thu, 30 Jul 2026, Rich Felker wrote:
> On Thu, Jul 30, 2026 at 10:16:49AM +0200, Paul Zimmermann wrote:
>> Hi,
>>
>> for x=-0x1.34e729fd08086p+21, Musl 1.2.6 yields 0x1.075fa39542c7fp+3
>> instead of NaN:
>>
>> $ VERBOSE=-v ./doit.musl acosh 1000
>> Checking acosh with musl-1.2.6
>> Using seed 16199
>> tested 1121438 numbers from extra file
>> Using 20 thread(s)
>> acosh 0 -1 -0x1.34e729fd08086p+21 [0] [inf] inf inf
>> libm gives 0x1.075fa39542c7fp+3
>> mpfr gives -nan
>
> Thanks! I looked at the function briefly. Does the comment:
>
> /* x < 1 domain error is handled in the called functions */
>
> somehow fail to hold because of a rounding error in the intermediate
> expression?
Yes.
I would explicitly handle the case of x < 1 to avoid such problems.
I can massage this to MUSL style if you think the fix makes sense.
double rcosh(double x)
{
static const unsigned long einf = 0x7ff;
static const unsigned long b = 0x3ff;
static const double negone = -1.0;
const union { double f; uint64_t i;} u = { x };
const uint64_t e = u.i >> 52; /* capture both sign bit plus exponent */
if (einf <= e || e < b ) /* x < 0 || x is a NaN || x == INF || x < 1 */
{
return x != x || x == INFINITY ? x + x : (x - x) / (x - x);
}
else if (e == b) /* 1 <= x < 2 */
{
return (x += negone, log1p(x + sqrt(2 * x + x * x)));
}
else if (e < b + 26) /* not too big */
{
return log((x + x) + negone/(x+sqrt(x * x + negone)));
}
return log(x) + 0.693147180559945309417232121458176568;
}
Regards - Damian
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
