This is the mail archive of the
libc-alpha@sources.redhat.com
mailing list for the glibc project.
nextafter and forced rounding
- From: Richard Henderson <rth at twiddle dot net>
- To: libc-alpha at sources dot redhat dot com
- Date: Fri, 6 Aug 2004 01:13:05 -0700
- Subject: nextafter and forced rounding
I was forced to look at nextafter implementations this week,
and I found we hadn't propagated the forced rounding fix to
all implementations.
The following pulls the messyness into a nice little inline
function and updates all of the implementations. I also add
some more nextafter test cases.
r~
2004-08-05 Richard Henderson <rth@redhat.com>
* math/libm-test.inc (nextafter_uf): New.
(nextafter_test): Test vs epsilon, and the minimum denormal.
(nexttoward_test): Likewise.
* math/math_private.h (__force_round_float): New.
(__force_round_double): New.
* sysdeps/generic/s_nextafter.c, sysdeps/generic/s_nexttowardf.c,
sysdeps/i386/fpu/s_nexttoward.c, sysdeps/i386/fpu/s_nexttowardf.c,
sysdeps/ieee754/flt-32/s_nextafterf.c,
sysdeps/ieee754/ldbl-128/s_nexttoward.c,
sysdeps/ieee754/ldbl-128/s_nexttowardf.c,
sysdeps/ieee754/ldbl-96/s_nexttoward.c,
sysdeps/ieee754/ldbl-96/s_nexttowardf.c: Use them.
Index: math/libm-test.inc
===================================================================
RCS file: /cvs/glibc/libc/math/libm-test.inc,v
retrieving revision 1.61
diff -u -p -r1.61 libm-test.inc
--- math/libm-test.inc 20 Jul 2004 07:06:32 -0000 1.61
+++ math/libm-test.inc 6 Aug 2004 07:52:15 -0000
@@ -3437,10 +3437,17 @@ nearbyint_test (void)
END (nearbyint);
}
+static FLOAT
+FUNC(nextafter_uf)(FLOAT x)
+{
+ x = x/2;
+ asm("" : "=m"(x) : "m"(x));
+ return x;
+}
+
static void
nextafter_test (void)
{
-
START (nextafter);
TEST_ff_f (nextafter, 0, 0, 0);
@@ -3462,6 +3469,16 @@ nextafter_test (void)
TEST_ff_f (nextafter, fltmax, plus_infty, plus_infty);
TEST_ff_f (nextafter, -fltmax, minus_infty, minus_infty);
+ FLOAT flteps = CHOOSE (LDBL_EPSILON, DBL_EPSILON, FLT_EPSILON,
+ LDBL_EPSILON, DBL_EPSILON, FLT_EPSILON);
+ TEST_ff_f (nextafter, 1, 2, 1+flteps);
+ TEST_ff_f (nextafter, 1+flteps, 0, 1);
+
+ FLOAT fltmin = FUNC(nextafter)(0, 1);
+ TEST_ff_f (nextafter, fltmin, -1, 0);
+ TEST_ff_f (nextafter, 0, -1, -fltmin);
+ TEST_f_f (nextafter_uf, fltmin, 0);
+
#ifdef TEST_LDOUBLE
// XXX Enable once gcc is fixed.
//TEST_ff_f (nextafter, 0x0.00000040000000000000p-16385L, -0.1L, 0x0.0000003ffffffff00000p-16385L);
@@ -3492,6 +3509,21 @@ nexttoward_test (void)
TEST_ff_f (nexttoward, 1.1L, nan_value, nan_value);
TEST_ff_f (nexttoward, nan_value, nan_value, nan_value);
+ FLOAT fltmax = CHOOSE (LDBL_MAX, DBL_MAX, FLT_MAX,
+ LDBL_MAX, DBL_MAX, FLT_MAX);
+ TEST_ff_f (nexttoward, fltmax, plus_infty, plus_infty);
+ TEST_ff_f (nexttoward, -fltmax, minus_infty, minus_infty);
+
+ FLOAT flteps = CHOOSE (LDBL_EPSILON, DBL_EPSILON, FLT_EPSILON,
+ LDBL_EPSILON, DBL_EPSILON, FLT_EPSILON);
+ TEST_ff_f (nexttoward, 1, 2, 1+flteps);
+ TEST_ff_f (nexttoward, 1+flteps, 0, 1);
+
+ FLOAT fltmin = FUNC(nexttoward)(0, 1);
+ TEST_ff_f (nexttoward, fltmin, -1, 0);
+ TEST_ff_f (nexttoward, 0, -1, -fltmin);
+ TEST_f_f (nextafter_uf, fltmin, 0);
+
/* XXX We need the hexadecimal FP number representation here for further
tests. */
Index: math/math_private.h
===================================================================
RCS file: /cvs/glibc/libc/math/math_private.h,v
retrieving revision 1.18
diff -u -p -r1.18 math_private.h
--- math/math_private.h 28 Aug 2003 00:10:15 -0000 1.18
+++ math/math_private.h 6 Aug 2004 07:52:15 -0000
@@ -150,6 +150,26 @@ do { \
/* Get long double macros from a separate header. */
#include <math_ldbl.h>
+/* Some helper functions to force rounding to the given precision. */
+
+#include <float.h>
+
+static inline float
+__force_round_float (float x)
+{
+ if (!__builtin_constant_p (x) && FLT_EVAL_METHOD != 0)
+ asm ("" : "=m"(x) : "m"(x));
+ return x;
+}
+
+static inline double
+__force_round_double (double x)
+{
+ if (!__builtin_constant_p (x) && FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
+ asm ("" : "=m"(x) : "m"(x));
+ return x;
+}
+
/* ieee style elementary functions */
extern double __ieee754_sqrt (double);
extern double __ieee754_acos (double);
Index: sysdeps/generic/s_nextafter.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/generic/s_nextafter.c,v
retrieving revision 1.6
diff -u -p -r1.6 s_nextafter.c
--- sysdeps/generic/s_nextafter.c 7 Dec 2003 20:53:42 -0000 1.6
+++ sysdeps/generic/s_nextafter.c 6 Aug 2004 07:52:16 -0000
@@ -50,7 +50,7 @@ static char rcsid[] = "$NetBSD: s_nextaf
if(x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
INSERT_WORDS(x,hy&0x80000000,1); /* return +-minsubnormal */
- y = x*x;
+ y = __force_round_double(x*x);
if(y==x) return y; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -71,14 +71,10 @@ static char rcsid[] = "$NetBSD: s_nextaf
}
}
hy = hx&0x7ff00000;
- if(hy>=0x7ff00000) {
- x = x+x; /* overflow */
- if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
- asm ("" : "=m"(x) : "m"(x));
- return x; /* overflow */
- }
+ if(hy>=0x7ff00000)
+ return __force_round_double(x+x); /* overflow */
if(hy<0x00100000) { /* underflow */
- y = x*x;
+ y = __force_round_double(x*x);
if(y!=x) { /* raise underflow flag */
INSERT_WORDS(y,hx,lx);
return y;
Index: sysdeps/generic/s_nexttowardf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/generic/s_nexttowardf.c,v
retrieving revision 1.5
diff -u -p -r1.5 s_nexttowardf.c
--- sysdeps/generic/s_nexttowardf.c 10 Dec 2003 22:52:08 -0000 1.5
+++ sysdeps/generic/s_nexttowardf.c 6 Aug 2004 07:52:16 -0000
@@ -47,7 +47,7 @@
if(ix==0) { /* x == 0 */
float x2;
SET_FLOAT_WORD(x,(u_int32_t)(hy&0x80000000)|1);/* return +-minsub*/
- x2 = x*x;
+ x2 = __force_round_float(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -66,15 +66,10 @@
hx += 1;
}
hy = hx&0x7f800000;
- if(hy>=0x7f800000) {
- x = x+x; /* overflow */
- if (FLT_EVAL_METHOD != 0)
- /* Force conversion to float. */
- asm ("" : "=m"(x) : "m"(x));
- return x;
- }
+ if(hy>=0x7f800000)
+ return __force_round_float(x+x); /* overflow */
if(hy<0x00800000) { /* underflow */
- float x2 = x*x;
+ float x2 = __force_round_float(x*x);
if(x2!=x) { /* raise underflow flag */
SET_FLOAT_WORD(x2,hx);
return x2;
Index: sysdeps/i386/fpu/s_nexttoward.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/i386/fpu/s_nexttoward.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_nexttoward.c
--- sysdeps/i386/fpu/s_nexttoward.c 7 Dec 2003 21:11:52 -0000 1.2
+++ sysdeps/i386/fpu/s_nexttoward.c 6 Aug 2004 07:52:16 -0000
@@ -54,7 +54,7 @@ static char rcsid[] = "$NetBSD: $";
if((ix|lx)==0) { /* x == 0 */
double x2;
INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
- x2 = x*x;
+ x2 = __force_round_double(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -83,14 +83,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7ff00000;
- if(hy>=0x7ff00000) {
- x = x+x; /* overflow */
- /* Force conversion to double. */
- asm ("" : "=m"(x) : "m"(x));
- return x;
- }
+ if(hy>=0x7ff00000)
+ return __force_round_double(x+x); /* overflow */
if(hy<0x00100000) { /* underflow */
- double x2 = x*x;
+ double x2 = __force_round_double(x*x);
if(x2!=x) { /* raise underflow flag */
INSERT_WORDS(x2,hx,lx);
return x2;
Index: sysdeps/i386/fpu/s_nexttowardf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/i386/fpu/s_nexttowardf.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_nexttowardf.c
--- sysdeps/i386/fpu/s_nexttowardf.c 7 Dec 2003 21:12:13 -0000 1.2
+++ sysdeps/i386/fpu/s_nexttowardf.c 6 Aug 2004 07:52:16 -0000
@@ -46,7 +46,7 @@ static char rcsid[] = "$NetBSD: $";
if(ix==0) { /* x == 0 */
float x2;
SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
- x2 = x*x;
+ x2 = __force_round_float(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -67,14 +67,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7f800000;
- if(hy>=0x7f800000) {
- x = x+x; /* overflow */
- /* Force conversion to float. */
- asm ("" : "=m"(x) : "m"(x));
- return x;
- }
+ if(hy>=0x7f800000)
+ return __force_round_float(x+x); /* overflow */
if(hy<0x00800000) { /* underflow */
- float x2 = x*x;
+ float x2 = __force_round_float(x*x);
if(x2!=x) { /* raise underflow flag */
SET_FLOAT_WORD(x2,hx);
return x2;
Index: sysdeps/ieee754/flt-32/s_nextafterf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/flt-32/s_nextafterf.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_nextafterf.c
--- sysdeps/ieee754/flt-32/s_nextafterf.c 7 Dec 2003 20:53:59 -0000 1.3
+++ sysdeps/ieee754/flt-32/s_nextafterf.c 6 Aug 2004 07:52:16 -0000
@@ -41,7 +41,7 @@ static char rcsid[] = "$NetBSD: s_nextaf
if(x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
- y = x*x;
+ y = __force_round_float(x*x);
if(y==x) return y; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -58,14 +58,10 @@ static char rcsid[] = "$NetBSD: s_nextaf
}
}
hy = hx&0x7f800000;
- if(hy>=0x7f800000) {
- x = x+x; /* overflow */
- if (FLT_EVAL_METHOD != 0)
- asm ("" : "=m"(x) : "m"(x));
- return x; /* overflow */
- }
+ if(hy>=0x7f800000)
+ return __force_round_float(x+x); /* overflow */
if(hy<0x00800000) { /* underflow */
- y = x*x;
+ y = __force_round_float(x*x);
if(y!=x) { /* raise underflow flag */
SET_FLOAT_WORD(y,hx);
return y;
Index: sysdeps/ieee754/ldbl-128/s_nexttoward.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-128/s_nexttoward.c,v
retrieving revision 1.3
diff -u -p -r1.3 s_nexttoward.c
--- sysdeps/ieee754/ldbl-128/s_nexttoward.c 7 Dec 2003 21:21:29 -0000 1.3
+++ sysdeps/ieee754/ldbl-128/s_nexttoward.c 6 Aug 2004 07:52:16 -0000
@@ -55,7 +55,7 @@ static char rcsid[] = "$NetBSD: $";
if((ix|lx)==0) { /* x == 0 */
double x2;
INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
- x2 = x*x;
+ x2 = __force_round_double(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -84,15 +84,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7ff00000;
- if(hy>=0x7ff00000) {
- x = x+x; /* overflow */
- if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
- /* Force conversion to float. */
- asm ("" : "=m"(x) : "m"(x));
- return x;
- }
+ if(hy>=0x7ff00000)
+ return __force_round_double(x+x); /* overflow */
if(hy<0x00100000) { /* underflow */
- double x2 = x*x;
+ double x2 = __force_round_double(x*x);
if(x2!=x) { /* raise underflow flag */
INSERT_WORDS(x2,hx,lx);
return x2;
Index: sysdeps/ieee754/ldbl-128/s_nexttowardf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-128/s_nexttowardf.c,v
retrieving revision 1.1
diff -u -p -r1.1 s_nexttowardf.c
--- sysdeps/ieee754/ldbl-128/s_nexttowardf.c 14 Jul 1999 00:09:42 -0000 1.1
+++ sysdeps/ieee754/ldbl-128/s_nexttowardf.c 6 Aug 2004 07:52:16 -0000
@@ -46,7 +46,7 @@ static char rcsid[] = "$NetBSD: $";
if(ix==0) { /* x == 0 */
float x2;
SET_FLOAT_WORD(x,(u_int32_t)((hy>>32)&0x80000000)|1);/* return +-minsub*/
- x2 = x*x;
+ x2 = __force_round_float(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -67,9 +67,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7f800000;
- if(hy>=0x7f800000) return x+x; /* overflow */
+ if(hy>=0x7f800000)
+ return __force_round_float(x+x); /* overflow */
if(hy<0x00800000) { /* underflow */
- float x2 = x*x;
+ float x2 = __force_round_float(x*x);
if(x2!=x) { /* raise underflow flag */
SET_FLOAT_WORD(x2,hx);
return x2;
Index: sysdeps/ieee754/ldbl-96/s_nexttoward.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-96/s_nexttoward.c,v
retrieving revision 1.5
diff -u -p -r1.5 s_nexttoward.c
--- sysdeps/ieee754/ldbl-96/s_nexttoward.c 7 Dec 2003 21:21:10 -0000 1.5
+++ sysdeps/ieee754/ldbl-96/s_nexttoward.c 6 Aug 2004 07:52:16 -0000
@@ -52,7 +52,7 @@ static char rcsid[] = "$NetBSD: $";
if((ix|lx)==0) { /* x == 0 */
double x2;
INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
- x2 = x*x;
+ x2 = __force_round_double(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -81,15 +81,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7ff00000;
- if(hy>=0x7ff00000) {
- x = x+x; /* overflow */
- if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
- /* Force conversion to float. */
- asm ("" : "=m"(x) : "m"(x));
- return x;
- }
+ if(hy>=0x7ff00000)
+ return __force_round_double(x+x); /* overflow */
if(hy<0x00100000) { /* underflow */
- double x2 = x*x;
+ double x2 = __force_round_double(x*x);
if(x2!=x) { /* raise underflow flag */
INSERT_WORDS(x2,hx,lx);
return x2;
Index: sysdeps/ieee754/ldbl-96/s_nexttowardf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/ldbl-96/s_nexttowardf.c,v
retrieving revision 1.2
diff -u -p -r1.2 s_nexttowardf.c
--- sysdeps/ieee754/ldbl-96/s_nexttowardf.c 16 Jun 2001 03:34:42 -0000 1.2
+++ sysdeps/ieee754/ldbl-96/s_nexttowardf.c 6 Aug 2004 07:52:16 -0000
@@ -43,7 +43,7 @@ static char rcsid[] = "$NetBSD: $";
if(ix==0) { /* x == 0 */
float x2;
SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
- x2 = x*x;
+ x2 = __force_round_float(x*x);
if(x2==x) return x2; else return x; /* raise underflow flag */
}
if(hx>=0) { /* x > 0 */
@@ -64,9 +64,10 @@ static char rcsid[] = "$NetBSD: $";
}
}
hy = hx&0x7f800000;
- if(hy>=0x7f800000) return x+x; /* overflow */
+ if(hy>=0x7f800000)
+ return __force_round_float (x+x); /* overflow */
if(hy<0x00800000) { /* underflow */
- float x2 = x*x;
+ float x2 = __force_round_float(x*x);
if(x2!=x) { /* raise underflow flag */
SET_FLOAT_WORD(x2,hx);
return x2;