Re: nextafterl(3) possible bug

[Daniel Dickman]

 confirming that this patch fixes the failing numpy regress test on i386.

 let me know if you want me to test a different diff.

 Here's a better diff, inspired by what FreeBSD has.

 ok?



ok with me. numpy works with this diff too.

Re: nextafterl(3) possible bug

[Mark Kettenis]

 Date: Mon, 2 Jun 2014 21:18:26 -0400
 From: Daniel Dickman didick...@gmail.com
 
 
  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.
 
 confirming that this patch fixes the failing numpy regress test on i386.
 
 let me know if you want me to test a different diff.

Here's a better diff, inspired by what FreeBSD has.

ok?


Index: 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
--- s_nextafterl.c  12 Nov 2013 21:07:28 -  1.4
+++ s_nextafterl.c  4 Jun 2014 10:05:17 -
@@ -32,8 +32,8 @@ nextafterl(long double x, long double y)
ix = esx0x7fff;/* |x| */
iy = esy0x7fff;/* |y| */
 
-   if (((ix==0x7fff)((hx|lx)!=0)) ||   /* x is nan */
-   ((iy==0x7fff)((hy|ly)!=0))) /* y is nan */
+   if (((ix==0x7fff)((hx0x7fff|lx)!=0)) ||   /* x is nan */
+   ((iy==0x7fff)((hy0x7fff|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,30 @@ nextafterl(long double x, long double y)
if(ixiy||((ix==iy)  (hxhy||((hx==hy)(lxly) {
  /* x  y, x -= ulp */
if(lx==0) {
-   if (hx==0) esx -= 1;
-   hx -= 1;
+   if ((hx0x7fff)==0) esx -= 1;
+   hx = (hx - 1) | (hx  0x8000);
}
lx -= 1;
} else {/* x  y, x += ulp */
lx += 1;
if(lx==0) {
-   hx += 1;
-   if (hx==0)
-   esx += 1;
+   hx = (hx + 1) | (hx  0x8000);
+   if ((hx0x7fff)==0) esx += 1;
}
}
} else {/* x  0 */
if(esy=0||(ixiy||((ix==iy)(hxhy||((hx==hy)(lxly)){
  /* x  y, x -= ulp */
if(lx==0) {
-   if (hx==0) esx -= 1;
-   hx -= 1;
+   if ((hx0x7fff)==0) esx -= 1;
+   hx = (hx - 1) | (hx  0x8000);
}
lx -= 1;
} else {/* x  y, x += ulp */
lx += 1;
if(lx==0) {
-   hx += 1;
-   if (hx==0) esx += 1;
+   hx = (hx + 1) | (hx  0x8000);
+   if ((hx0x7fff)==0) esx += 1;
}
}
}

nextafterl(3) possible bug

[Daniel Dickman]

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.

---[ on linux/x86_64 ]---

# gcc -lm test2.c  ./a.out
1.00e+00
-5.421011e-20
OK
1.00
-0.00
OK

---[ on openbsd/i386 ]---

# gcc -lm test2.c  ./a.out
1.00e+00
nan   -- this looks wrong
1.00
-0.00
OK

---[ test2.c ]---

#include stdio.h
#include stdlib.h
#include math.h

int main(void) {
long double one_ld;
long double zero_ld;
long double after_ld;
long double result_ld;

one_ld = 1.0L;
zero_ld = 0.0L;

after_ld = nextafterl(one_ld, zero_ld);
printf(%Le\n, after_ld);
result_ld = after_ld - one_ld;
// result_ld should not be nan
printf(%Le\n, result_ld);
if (result_ld  0)
printf(OK\n);

float one_f;
float zero_f;
float after_f;
float result_f;

one_f = 1.0;
zero_f = 0.0;

after_f = nextafterf(one_f, zero_f);
printf(%f\n, after_f);
result_f = after_f - one_f;
printf(%f\n, result_f);
if (result_f  0)
printf(OK\n);
}

Re: nextafterl(3) possible bug

[Mark Kettenis]

 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 -  1.4
+++ src/ld80/s_nextafterl.c 2 Jun 2014 13:21:58 -
@@ -32,8 +32,8 @@ nextafterl(long double x, long double y)
ix = esx0x7fff;/* |x| */
iy = esy0x7fff;/* |y| */
 
-   if (((ix==0x7fff)((hx|lx)!=0)) ||   /* x is nan */
-   ((iy==0x7fff)((hy|ly)!=0))) /* y is nan */
+   if (((ix==0x7fff)((hx=0x7fff|lx)!=0)) ||   /* x is nan */
+   ((iy==0x7fff)((hy-0x7fff|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(ixiy||((ix==iy)  (hxhy||((hx==hy)(lxly) {
  /* x  y, x -= ulp */
if(lx==0) {
-   if (hx==0) esx -= 1;
-   hx -= 1;
+   if (hx = 0x8000) {
+ if (esx == 0) {
+   --hx;
+ } else {
+   esx -= 1;
+   hx = hx - 1;
+   if (esx  0)
+ hx |= 0x8000;
+ }
+   } 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 == 0x8000)) {
esx += 1;
+   hx |= 0x8000;
+   }
}
}
} else {/* x  0 */
if(esy=0||(ixiy||((ix==iy)(hxhy||((hx==hy)(lxly)){
  /* x  y, x -= ulp */
if(lx==0) {
-   if (hx==0) esx -= 1;
-   hx -= 1;
+   if (hx = 0x8000) {
+   esx -= 1;
+   hx = hx - 1;
+   if ((esx0x7fff)  0)
+ hx |= 0x8000;
+   } 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 == 0x8000  hx == 0x8000)) {
+   esx += 1;
+   hx |= 0x8000;
+   }
}
}
}

Re: nextafterl(3) possible bug

[Daniel Dickman]

 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.

confirming that this patch fixes the failing numpy regress test on i386.

let me know if you want me to test a different diff.


 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 -  1.4
 +++ src/ld80/s_nextafterl.c 2 Jun 2014 13:21:58 -
 @@ -32,8 +32,8 @@ nextafterl(long double x, long double y)
 ix = esx0x7fff;/* |x| */
 iy = esy0x7fff;/* |y| */

 -   if (((ix==0x7fff)((hx|lx)!=0)) ||   /* x is nan */
 -   ((iy==0x7fff)((hy|ly)!=0))) /* y is nan */
 +   if (((ix==0x7fff)((hx=0x7fff|lx)!=0)) ||   /* x is nan */
 +   ((iy==0x7fff)((hy-0x7fff|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(ixiy||((ix==iy)  (hxhy||((hx==hy)(lxly) {
   /* x  y, x -= ulp */
 if(lx==0) {
 -   if (hx==0) esx -= 1;
 -   hx -= 1;
 +   if (hx = 0x8000) {
 + if (esx == 0) {
 +   --hx;
 + } else {
 +   esx -= 1;
 +   hx = hx - 1;
 +   if (esx  0)
 + hx |= 0x8000;
 + }
 +   } 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 == 0x8000)) {
 esx += 1;
 +   hx |= 0x8000;
 +   }
 }
 }
 } else {/* x  0 */
 if(esy=0||(ixiy||((ix==iy)(hxhy||((hx==hy)(lxly)){
   /* x  y, x -= ulp */
 if(lx==0) {
 -   if (hx==0) esx -= 1;
 -   hx -= 1;
 +   if (hx = 0x8000) {
 +   esx -= 1;
 +   hx = hx - 1;
 +   if ((esx0x7fff)  0)
 + hx |= 0x8000;
 +   } 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 == 0x8000  hx == 0x8000)) {
 +   esx += 1;
 +   hx |= 0x8000;
 +   }
 }
 }
 }