> Date: Mon, 2 Jun 2014 07:34:59 -0400
> From: Daniel Dickman <didick...@gmail.com>
>
> >From the numpy test suite, I think I might have found a bug in
> nextafterl(3). The "result_ld" variable below comes back as nan on
> i386. But doing the same calculations with floats returns the expected
> values.
>
> A test on Linux also shows the expected results for both the float and
> long double cases.
Another bug. Intel chose an extended precision format with an
explicit integer bit, and the code doesn't handle that. Assuming we
don't support machines with extended precision format that have an
implicit integer bit, the following diff (an adaptation of the code in
glibc) should fix things. Not entirely happy with the fix though, so
I'm still thinking about improvements.
Index: src/ld80/s_nextafterl.c
===================================================================
RCS file: /cvs/src/lib/libm/src/ld80/s_nextafterl.c,v
retrieving revision 1.4
diff -u -p -r1.4 s_nextafterl.c
--- src/ld80/s_nextafterl.c 12 Nov 2013 21:07:28 -0000 1.4
+++ src/ld80/s_nextafterl.c 2 Jun 2014 13:21:58 -0000
@@ -32,8 +32,8 @@ nextafterl(long double x, long double y)
ix = esx&0x7fff; /* |x| */
iy = esy&0x7fff; /* |y| */
- if (((ix==0x7fff)&&((hx|lx)!=0)) || /* x is nan */
- ((iy==0x7fff)&&((hy|ly)!=0))) /* y is nan */
+ if (((ix==0x7fff)&&((hx&=0x7fffffff|lx)!=0)) || /* x is nan */
+ ((iy==0x7fff)&&((hy&-0x7fffffff|ly)!=0))) /* y is nan */
return x+y;
if(x==y) return y; /* x=y, return y */
if((ix|hx|lx)==0) { /* x == 0 */
@@ -47,31 +47,50 @@ nextafterl(long double x, long double y)
if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
/* x > y, x -= ulp */
if(lx==0) {
- if (hx==0) esx -= 1;
- hx -= 1;
+ if (hx <= 0x80000000) {
+ if (esx == 0) {
+ --hx;
+ } else {
+ esx -= 1;
+ hx = hx - 1;
+ if (esx > 0)
+ hx |= 0x80000000;
+ }
+ } else
+ hx -= 1;
}
lx -= 1;
} else { /* x < y, x += ulp */
lx += 1;
if(lx==0) {
hx += 1;
- if (hx==0)
+ if (hx==0 || (esx == 0 && hx == 0x80000000)) {
esx += 1;
+ hx |= 0x80000000;
+ }
}
}
} else { /* x < 0 */
if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){
/* x < y, x -= ulp */
if(lx==0) {
- if (hx==0) esx -= 1;
- hx -= 1;
+ if (hx <= 0x80000000) {
+ esx -= 1;
+ hx = hx - 1;
+ if ((esx&0x7fff) > 0)
+ hx |= 0x80000000;
+ } else
+ hx -= 1;
}
lx -= 1;
} else { /* x > y, x += ulp */
lx += 1;
if(lx==0) {
hx += 1;
- if (hx==0) esx += 1;
+ if (hx==0 || (esx == 0xffff8000 && hx == 0x80000000)) {
+ esx += 1;
+ hx |= 0x80000000;
+ }
}
}
}
