[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <20260731005938.GK27423@brightrain.aerifal.cx>
Date: Thu, 30 Jul 2026 20:59:38 -0400
From: Rich Felker <dalias@...c.org>
To: Damian McGuckin <damianm@....com.au>
Cc: musl@...ts.openwall.com, Paul Zimmermann <Paul.Zimmermann@...ia.fr>
Subject: Re: issue in acosh
On Fri, Jul 31, 2026 at 10:01:38AM +1000, Damian McGuckin wrote:
> 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);
> }
I don't think nan/inf need special treatment here. Just
if (e < 0x3ff) x = 0;
should make it work fine.
Rich
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.
