Follow @Openwall on Twitter for new release announcements and other news
[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Date: Sun, 15 Aug 2021 10:03:12 -0400
From: Rich Felker <dalias@...c.org>
To: Damian McGuckin <damianm@....com.au>
Cc: musl@...ts.openwall.com, Szabolcs Nagy <nsz@...t70.net>
Subject: Re: [PATCH #2] Properly simplified nextafter()

On Sun, Aug 15, 2021 at 06:24:31PM +1000, Damian McGuckin wrote:
> 
> Hi Stefan,
> 
> On Sun, 15 Aug 2021, Stefan Kanthak wrote:
> 
> >__attribute__((noinline))
> >double nextafter(double x, double y)
> >{
> >union {double f; unsigned long long i;} ux={x}, uy={y};
> >unsigned long long ax, ay;
> >int e;
> >
> >if (isnan(x) || isnan(y))
> > return x + y;
> >if (ux.i == uy.i)
> > return y;
> >#ifdef PATCH
> >ax = ux.i << 1;
> >ay = uy.i << 1;
> >#else
> >ax = ux.i & -1ULL/2;
> >ay = uy.i & -1ULL/2;
> >#endif
> >if (ax == 0) {
> > if (ay == 0)
> >  return y;
> > ux.i = (uy.i & 1ULL<<63) | 1;
> >#ifdef PATCH
> >} else if (ax < ay == (long long) ux.i < 0)
> >#else
> >} else if (ax > ay || ((ux.i ^ uy.i) & 1ULL<<63))
> >#endif
> > ux.i--;
> >else
> > ux.i++;
> >e = ux.i >> 52 & 0x7ff;
> >/* raise overflow if ux.f is infinite and x is finite */
> >if (e == 0x7ff)
> > FORCE_EVAL(x + x);
> >/* raise underflow if ux.f is subnormal or zero */
> >if (e == 0)
> > FORCE_EVAL(x*x + ux.f*ux.f);
> >return ux.f;
> >}
> 
> Maybe I am missing something and my brain is in weekend-mode ...
> 
> I did a quick check and ran the above code for some test cases:
> 
> nextafter(-9.7500000000e+01, 3.5000000000e+01) = -9.7500000000e+01
> yourpatch(-9.7500000000e+01, 3.5000000000e+01) = -9.7500000000e+01
> 
> The error is 2.8421709430e-14
> 
> nextafter(-9.7500000000e+01, -3.5000000000e+01) = -9.7500000000e+01
> yourpatch(-9.7500000000e+01, -3.5000000000e+01) = -9.7500000000e+01
> 
> The error is 2.8421709430e-14

I don't follow; are you claiming Stefan's patch introduces an error
here? The outputs you printed show the exact same behavior before and
after but it's possible you printed them wrong. The nextafter function
is bit-exact; it does not have floating point error (inexactness)
unless the implementation is buggy.

> nextafter(-inf, inf) = -1.7976931349e+308 Correct
> yourpatch(-inf, inf) = -nan Incorrect
> 
> This is against standard GLIB.

glibc I assume you mean? In any case yes this looks like a bug in the
patch.

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.