Re: [musl] [PATCH #2] Properly simplified nextafter()

[Rich Felker <dalias@libc.org>]

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