public inbox for libc-hacker@sourceware.org
 help / color / mirror / Atom feed
From: Jakub Jelinek <jakub@redhat.com>
To: Ulrich Drepper <drepper@redhat.com>
Cc: Glibc hackers <libc-hacker@sources.redhat.com>
Subject: [PATCH] Fix next{after,toward}{,f,l} underflow (and overflow) handling
Date: Tue, 27 Mar 2007 13:08:00 -0000	[thread overview]
Message-ID: <20070327131242.GD1826@sunsite.mff.cuni.cz> (raw)

Hi!

Apparently GCC will DCE even floating point operations that might raise
exceptions and clearly has been doing that for years.
But POSIX requires exceptions to be raised in certain situations
from nextafter*/nexttoward*, particularly if x != y and the return value is 0
or subnormal, FE_UNDERFLOW should be raised and if FE_OVERFLOW if it is
returning +-Inf.

The generic math_opt_barrier and math_force_eval definitions ought to work
everywhere, but on some arches might pessimize the code.  The first macro
should just return its argument and hide it from the optimizer, so generally
doesn't need to go through mem, if some arch has special constraints for
floating point regs, it should use them.  There is no need for this
macro to force rounding to the actual type.  math_force_eval needs to ensure
rounding to the actual type of the var and then the result can be thrown
away, all we care about is side-effects.

2007-03-27  Jakub Jelinek  <jakub@redhat.com>

	[BZ #3306]
	* math/math_private.h (math_opt_barrier, math_force_eval): Define.
	* sysdeps/i386/fpu/math_private.h: New file.
	* sysdeps/x86_64/fpu/math_private.h: New file.
	* math/s_nexttowardf.c (__nexttowardf): Use math_opt_barrier and
	math_force_eval macros.  Use "+m" constraint on asm rather than
	"=m" and "m".
	* math/s_nextafter.c (__nextafter): Likewise.
	* sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c (__nexttoward):
	Likewise.
	* sysdeps/ieee754/flt-32/s_nextafterf.c (__nextafterf): Likewise.
	* sysdeps/ieee754/ldbl-128/s_nexttoward.c (__nexttoward): Likewise.
	* sysdeps/ieee754/ldbl-96/s_nexttoward.c (__nexttoward): Likewise.
	* sysdeps/i386/fpu/s_nextafterl.c (__nextafterl): Use
	math_opt_barrier and math_force_eval macros.
	* sysdeps/ieee754/ldbl-128/s_nextafterl.c (__nextafterl): Likewise.
	* sysdeps/ieee754/ldbl-96/s_nextafterl.c (__nextafterl): Likewise.
	* sysdeps/i386/fpu/s_nexttoward.c: Include float.h.
	(__nexttoward): Use math_opt_barrier and
	math_force_eval macros.  Use "+m" constraint on asm rather than
	"=m" and "m".  Only use asm to force double result if
	FLT_EVAL_METHOD is 2.
	* sysdeps/i386/fpu/s_nexttowardf.c: Include float.h.
	(__nexttowardf): Use math_opt_barrier and
	math_force_eval macros.  Use "+m" constraint on asm rather than
	"=m" and "m".  Only use asm to force double result if
	FLT_EVAL_METHOD is not 0.
	* sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c: Include float.h.
	(__nexttowardf): Use math_opt_barrier and
	math_force_eval macros.  If FLT_EVAL_METHOD is not 0, force
	x to float using asm.
	* sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c: Include float.h.
	(__nldbl_nexttowardf): Use math_opt_barrier and
	math_force_eval macros.  If FLT_EVAL_METHOD is not 0, force
	x to float using asm.
	* sysdeps/ieee754/ldbl-96/s_nexttowardf.c: Include float.h.
	(__nexttowardf): Use math_opt_barrier and math_force_eval
	macros.  If FLT_EVAL_METHOD is not 0, force x to float using asm.
	* math/bug-nextafter.c (zero, inf): New variables.
	(main): Add new tests.
	* math/bug-nexttoward.c (zero, inf): New variables.
	(main): Add new tests.

--- libc/math/s_nexttowardf.c.jj	2005-12-14 11:34:14.000000000 +0100
+++ libc/math/s_nexttowardf.c	2007-03-27 14:05:13.000000000 +0200
@@ -21,7 +21,7 @@
  */
 
 #include <math.h>
-#include "math_private.h"
+#include <math_private.h>
 #include <float.h>
 
 #ifdef __STDC__
@@ -45,10 +45,12 @@
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
-	    float x2;
+	    float u;
 	    SET_FLOAT_WORD(x,(u_int32_t)(hy&0x80000000)|1);/* return +-minsub*/
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x; /* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		 /* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hy<0||(ix>>23)>(iy>>20)-0x380
@@ -70,15 +72,12 @@
 	  x = x+x;	/* overflow  */
 	  if (FLT_EVAL_METHOD != 0)
 	    /* Force conversion to float.  */
-	    asm ("" : "=m"(x) : "m"(x));
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00800000) {		/* underflow */
-	    float x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-		SET_FLOAT_WORD(x2,hx);
-		return x2;
-	    }
+	if(hy<0x00800000) {
+	    float u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/math/bug-nextafter.c.jj	2003-12-09 19:10:47.000000000 +0100
+++ libc/math/bug-nextafter.c	2007-03-23 15:25:00.000000000 +0100
@@ -4,6 +4,9 @@
 #include <stdlib.h>
 #include <stdio.h>
 
+float zero = 0.0;
+float inf = INFINITY;
+
 int
 main (void)
 {
@@ -34,6 +37,81 @@ main (void)
       ++result;
     }
 
+  i = 0;
+  m = FLT_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  i = nextafterf (m, i);
+  if (i < 0 || i >= FLT_MIN)
+    {
+      puts ("nextafterf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf+ did not underflow");
+      ++result;
+    }
+  i = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  i = nextafterf (-m, -i);
+  if (i > 0 || i <= -FLT_MIN)
+    {
+      puts ("nextafterf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf- did not underflow");
+      ++result;
+    }
+  i = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  m = nextafterf (zero, inf);
+  if (m < 0.0 || m >= FLT_MIN)
+    {
+      puts ("nextafterf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterf (m, i) != 0.0)
+    {
+      puts ("nextafterf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  m = nextafterf (copysignf (zero, -1.0), -inf);
+  if (m > 0.0 || m <= -FLT_MIN)
+    {
+      puts ("nextafterf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterf (m, -i) != 0.0)
+    {
+      puts ("nextafterf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterf- did not underflow");
+      ++result;
+    }
+
   double di = INFINITY;
   double dm = DBL_MAX;
   feclearexcept (FE_ALL_EXCEPT);
@@ -59,5 +137,182 @@ main (void)
       ++result;
     }
 
+  di = 0;
+  dm = DBL_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  di = nextafter (dm, di);
+  if (di < 0 || di >= DBL_MIN)
+    {
+      puts ("nextafter+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter+ did not underflow");
+      ++result;
+    }
+  di = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  di = nextafter (-dm, -di);
+  if (di > 0 || di <= -DBL_MIN)
+    {
+      puts ("nextafter- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter- did not underflow");
+      ++result;
+    }
+  di = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  dm = nextafter (zero, inf);
+  if (dm < 0.0 || dm >= DBL_MIN)
+    {
+      puts ("nextafter+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafter (dm, di) != 0.0)
+    {
+      puts ("nextafter+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  dm = nextafter (copysign (zero, -1.0), -inf);
+  if (dm > 0.0 || dm <= -DBL_MIN)
+    {
+      puts ("nextafter- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafter (dm, -di) != 0.0)
+    {
+      puts ("nextafter- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafter- did not underflow");
+      ++result;
+    }
+
+#ifndef NO_LONG_DOUBLE
+  long double li = INFINITY;
+  long double lm = LDBL_MAX;
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterl (lm, li) != li)
+    {
+      puts ("nextafterl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_OVERFLOW) == 0)
+    {
+      puts ("nextafterl+ did not overflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterl (-lm, -li) != -li)
+    {
+      puts ("nextafterl failed");
+      ++result;
+    }
+  if (fetestexcept (FE_OVERFLOW) == 0)
+    {
+      puts ("nextafterl- did not overflow");
+      ++result;
+    }
+
+  li = 0;
+  lm = LDBL_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  li = nextafterl (lm, li);
+  if (li < 0 || li >= LDBL_MIN)
+    {
+      puts ("nextafterl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl+ did not underflow");
+      ++result;
+    }
+  li = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  li = nextafterl (-lm, -li);
+  if (li > 0 || li <= -LDBL_MIN)
+    {
+      puts ("nextafterl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl- did not underflow");
+      ++result;
+    }
+  li = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  lm = nextafterl (zero, inf);
+  if (lm < 0.0 || lm >= LDBL_MIN)
+    {
+      puts ("nextafterl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterl (lm, li) != 0.0)
+    {
+      puts ("nextafterl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  lm = nextafterl (copysign (zero, -1.0), -inf);
+  if (lm > 0.0 || lm <= -LDBL_MIN)
+    {
+      puts ("nextafterl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nextafterl (lm, -li) != 0.0)
+    {
+      puts ("nextafterl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nextafterl- did not underflow");
+      ++result;
+    }
+#endif
+
   return result;
 }
--- libc/math/s_nextafter.c.jj	2005-12-14 11:33:59.000000000 +0100
+++ libc/math/s_nextafter.c	2007-03-27 14:05:23.000000000 +0200
@@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: s_nextaf
 #define nexttoward __internal_nexttoward
 
 #include <math.h>
-#include "math_private.h"
+#include <math_private.h>
 #include <float.h>
 
 #ifdef __STDC__
@@ -49,9 +49,12 @@ static char rcsid[] = "$NetBSD: s_nextaf
 	   return x+y;
 	if(x==y) return y;		/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
+	    double u;
 	    INSERT_WORDS(x,hy&0x80000000,1);	/* return +-minsubnormal */
-	    y = x*x;
-	    if(y==x) return y; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u*u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
@@ -74,15 +77,12 @@ static char rcsid[] = "$NetBSD: s_nextaf
 	if(hy>=0x7ff00000) {
 	  x = x+x;	/* overflow  */
 	  if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
-	    asm ("" : "=m"(x) : "m"(x));
+	    asm ("" : "+m"(x));
 	  return x;	/* overflow  */
 	}
-	if(hy<0x00100000) {		/* underflow */
-	    y = x*x;
-	    if(y!=x) {		/* raise underflow flag */
-	        INSERT_WORDS(y,hx,lx);
-		return y;
-	    }
+	if(hy<0x00100000) {
+	    double u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	INSERT_WORDS(x,hx,lx);
 	return x;
--- libc/math/math_private.h.jj	2005-11-21 16:43:03.000000000 +0100
+++ libc/math/math_private.h	2007-03-23 14:32:06.000000000 +0100
@@ -332,4 +332,10 @@ extern double __slowexp (double __x);
 extern double __slowpow (double __x, double __y, double __z);
 extern void __docos (double __x, double __dx, double __v[]);
 
+#ifndef math_opt_barrier
+#define math_opt_barrier(x) \
+({ __typeof (x) __x = x; __asm ("" : "+m" (__x)); __x; })
+#define math_force_eval(x) __asm __volatile ("" : : "m" (x))
+#endif
+
 #endif /* _MATH_PRIVATE_H_ */
--- libc/math/bug-nexttoward.c.jj	2003-12-07 22:13:09.000000000 +0100
+++ libc/math/bug-nexttoward.c	2007-03-23 15:47:10.000000000 +0100
@@ -4,6 +4,9 @@
 #include <stdlib.h>
 #include <stdio.h>
 
+float zero = 0.0;
+float inf = INFINITY;
+
 int
 main (void)
 {
@@ -35,6 +38,81 @@ main (void)
       ++result;
     }
 
+  fi = 0;
+  m = FLT_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  fi = nexttowardf (m, fi);
+  if (fi < 0 || fi >= FLT_MIN)
+    {
+      puts ("nexttowardf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf+ did not underflow");
+      ++result;
+    }
+  fi = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  fi = nexttowardf (-m, -fi);
+  if (fi > 0 || fi <= -FLT_MIN)
+    {
+      puts ("nexttowardf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf- did not underflow");
+      ++result;
+    }
+  fi = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  m = nexttowardf (zero, inf);
+  if (m < 0.0 || m >= FLT_MIN)
+    {
+      puts ("nexttowardf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardf (m, fi) != 0.0)
+    {
+      puts ("nexttowardf+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  m = nexttowardf (copysignf (zero, -1.0), -inf);
+  if (m > 0.0 || m <= -FLT_MIN)
+    {
+      puts ("nexttowardf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardf (m, -fi) != 0.0)
+    {
+      puts ("nexttowardf- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardf- did not underflow");
+      ++result;
+    }
+
   tl = (long double) DBL_MAX + 1.0e305L;
   double di = INFINITY;
   double dm = DBL_MAX;
@@ -61,5 +139,182 @@ main (void)
       ++result;
     }
 
+  di = 0;
+  dm = DBL_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  di = nexttoward (dm, di);
+  if (di < 0 || di >= DBL_MIN)
+    {
+      puts ("nexttoward+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward+ did not underflow");
+      ++result;
+    }
+  di = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  di = nexttoward (-dm, -di);
+  if (di > 0 || di <= -DBL_MIN)
+    {
+      puts ("nexttoward- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward- did not underflow");
+      ++result;
+    }
+  di = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  dm = nexttoward (zero, inf);
+  if (dm < 0.0 || dm >= DBL_MIN)
+    {
+      puts ("nexttoward+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttoward (dm, di) != 0.0)
+    {
+      puts ("nexttoward+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  dm = nexttoward (copysign (zero, -1.0), -inf);
+  if (dm > 0.0 || dm <= -DBL_MIN)
+    {
+      puts ("nexttoward- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttoward (dm, -di) != 0.0)
+    {
+      puts ("nexttoward- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttoward- did not underflow");
+      ++result;
+    }
+
+#ifndef NO_LONG_DOUBLE
+  long double li = INFINITY;
+  long double lm = LDBL_MAX;
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardl (lm, li) != li)
+    {
+      puts ("nexttowardl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_OVERFLOW) == 0)
+    {
+      puts ("nexttowardl+ did not overflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardl (-lm, -li) != -li)
+    {
+      puts ("nexttowardl failed");
+      ++result;
+    }
+  if (fetestexcept (FE_OVERFLOW) == 0)
+    {
+      puts ("nexttowardl- did not overflow");
+      ++result;
+    }
+
+  li = 0;
+  lm = LDBL_MIN;
+  feclearexcept (FE_ALL_EXCEPT);
+  li = nexttowardl (lm, li);
+  if (li < 0 || li >= LDBL_MIN)
+    {
+      puts ("nexttowardl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl+ did not underflow");
+      ++result;
+    }
+  li = 0;
+  feclearexcept (FE_ALL_EXCEPT);
+  li = nexttowardl (-lm, -li);
+  if (li > 0 || li <= -LDBL_MIN)
+    {
+      puts ("nexttowardl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl- did not underflow");
+      ++result;
+    }
+  li = -INFINITY;
+  feclearexcept (FE_ALL_EXCEPT);
+  lm = nexttowardl (zero, inf);
+  if (lm < 0.0 || lm >= LDBL_MIN)
+    {
+      puts ("nexttowardl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardl (lm, li) != 0.0)
+    {
+      puts ("nexttowardl+ failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl+ did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  lm = nexttowardl (copysign (zero, -1.0), -inf);
+  if (lm > 0.0 || lm <= -LDBL_MIN)
+    {
+      puts ("nexttowardl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl- did not underflow");
+      ++result;
+    }
+  feclearexcept (FE_ALL_EXCEPT);
+  if (nexttowardl (lm, -li) != 0.0)
+    {
+      puts ("nexttowardl- failed");
+      ++result;
+    }
+  if (fetestexcept (FE_UNDERFLOW) == 0)
+    {
+      puts ("nexttowardl- did not underflow");
+      ++result;
+    }
+#endif
+
   return result;
 }
--- libc/sysdeps/i386/fpu/s_nextafterl.c.jj	2003-02-15 06:22:58.000000000 +0100
+++ libc/sysdeps/i386/fpu/s_nextafterl.c	2007-03-27 13:50:40.000000000 +0200
@@ -27,7 +27,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 
 #ifdef __STDC__
 	long double __nextafterl(long double x, long double y)
@@ -52,9 +52,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if(x==y) return y;		/* x=y, return y */
 	if((ix|hx|lx)==0) {			/* x == 0 */
+	    long double u;
 	    SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */
-	    y = x*x;
-	    if(y==x) return y; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(esx>=0) {			/* x > 0 */
 	    if(esx>esy||((esx==esy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
@@ -109,12 +112,9 @@ static char rcsid[] = "$NetBSD: $";
 	}
 	esy = esx&0x7fff;
 	if(esy==0x7fff) return x+x;	/* overflow  */
-	if(esy==0) {			/* underflow */
-	    y = x*x;
-	    if(y!=x) {		/* raise underflow flag */
-	        SET_LDOUBLE_WORDS(y,esx,hx,lx);
-		return y;
-	    }
+	if(esy==0) {
+	    long double u = x*x;		/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_LDOUBLE_WORDS(x,esx,hx,lx);
 	return x;
--- libc/sysdeps/i386/fpu/s_nexttoward.c.jj	2007-03-27 13:28:32.000000000 +0200
+++ libc/sysdeps/i386/fpu/s_nexttoward.c	2007-03-27 14:09:50.000000000 +0200
@@ -27,7 +27,8 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
+#include <float.h>
 
 #ifdef __STDC__
 	double __nexttoward(double x, long double y)
@@ -52,10 +53,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
-	    double x2;
+	    double u;
 	    INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00
@@ -85,16 +88,14 @@ static char rcsid[] = "$NetBSD: $";
 	hy = hx&0x7ff00000;
 	if(hy>=0x7ff00000) {
 	  x = x+x;	/* overflow  */
-	  /* Force conversion to double.  */
-	  asm ("" : "=m"(x) : "m"(x));
+	  if (FLT_EVAL_METHOD != 0 && FLT_EVAL_METHOD != 1)
+	    /* Force conversion to double.  */
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00100000) {		/* underflow */
-	    double x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        INSERT_WORDS(x2,hx,lx);
-		return x2;
-	    }
+	if(hy<0x00100000) {
+	    double u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	INSERT_WORDS(x,hx,lx);
 	return x;
--- libc/sysdeps/i386/fpu/s_nexttowardf.c.jj	2003-12-07 22:12:13.000000000 +0100
+++ libc/sysdeps/i386/fpu/s_nexttowardf.c	2007-03-27 14:09:58.000000000 +0200
@@ -19,7 +19,8 @@ static char rcsid[] = "$NetBSD: $";
 #endif
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
+#include <float.h>
 
 #ifdef __STDC__
 	float __nexttowardf(float x, long double y)
@@ -44,10 +45,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
-	    float x2;
+	    float u;
 	    SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80
@@ -69,16 +72,14 @@ static char rcsid[] = "$NetBSD: $";
 	hy = hx&0x7f800000;
 	if(hy>=0x7f800000) {
 	  x = x+x;	/* overflow  */
-	  /* Force conversion to float.  */
-	  asm ("" : "=m"(x) : "m"(x));
+	  if (FLT_EVAL_METHOD != 0)
+	    /* Force conversion to float.  */
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00800000) {		/* underflow */
-	    float x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        SET_FLOAT_WORD(x2,hx);
-		return x2;
-	    }
+	if(hy<0x00800000) {
+	    float u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/sysdeps/i386/fpu/math_private.h.jj	2007-03-23 14:30:22.000000000 +0100
+++ libc/sysdeps/i386/fpu/math_private.h	2007-03-23 15:32:38.000000000 +0100
@@ -0,0 +1,18 @@
+#ifndef _MATH_PRIVATE_H
+
+#define math_opt_barrier(x) \
+({ __typeof(x) __x;					\
+   __asm ("" : "=t" (__x) : "0" (x));			\
+   __x; })
+#define math_force_eval(x) \
+do							\
+  {							\
+    if (sizeof (x) <= sizeof (double))			\
+      __asm __volatile ("" : : "m" (x));		\
+    else						\
+      __asm __volatile ("" : : "f" (x));		\
+  }							\
+while (0)
+
+#include <math/math_private.h>
+#endif
--- libc/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c.jj	2006-01-28 01:07:25.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_nexttoward.c	2007-03-27 14:06:04.000000000 +0200
@@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 #include <math_ldbl_opt.h>
 #include <float.h>
 
@@ -55,10 +55,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
-	    double x2;
+	    double u;
 	    INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if (hy<0||(ix>>20)>(iy>>52)
@@ -89,16 +91,13 @@ static char rcsid[] = "$NetBSD: $";
 	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));
+	    /* Force conversion to double.  */
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00100000) {		/* underflow */
-	    double x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        INSERT_WORDS(x2,hx,lx);
-		return x2;
-	    }
+	if(hy<0x00100000) {
+	    double u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	INSERT_WORDS(x,hx,lx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c.jj	2006-01-28 01:07:25.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128ibm/s_nexttowardf.c	2007-03-27 14:11:11.000000000 +0200
@@ -19,8 +19,9 @@ static char rcsid[] = "$NetBSD: $";
 #endif
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 #include <math_ldbl_opt.h>
+#include <float.h>
 
 #ifdef __STDC__
 	float __nexttowardf(float x, long double y)
@@ -46,10 +47,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
-	    float x2;
+	    float u;
 	    SET_FLOAT_WORD(x,(u_int32_t)((hy>>32)&0x80000000)|1);/* return +-minsub*/
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hy<0||(ix>>23)>(iy>>52)-0x380
@@ -69,13 +72,16 @@ static char rcsid[] = "$NetBSD: $";
 	    }
 	}
 	hy = hx&0x7f800000;
-	if(hy>=0x7f800000) return x+x;	/* overflow  */
+	if(hy>=0x7f800000) {
+	  x = x+x;	/* overflow  */
+	  if (FLT_EVAL_METHOD != 0)
+	    /* Force conversion to float.  */
+	    asm ("" : "+m"(x));
+	  return x;
+	}
 	if(hy<0x00800000) {		/* underflow */
-	    float x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        SET_FLOAT_WORD(x2,hx);
-		return x2;
-	    }
+	    float u = x*x;
+	    math_force_eval (u);	/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/sysdeps/ieee754/flt-32/s_nextafterf.c.jj	2003-12-07 21:53:59.000000000 +0100
+++ libc/sysdeps/ieee754/flt-32/s_nextafterf.c	2007-03-27 14:06:23.000000000 +0200
@@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: s_nextaf
 #endif
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 #include <float.h>
 
 #ifdef __STDC__
@@ -40,9 +40,12 @@ static char rcsid[] = "$NetBSD: s_nextaf
 	   return x+y;
 	if(x==y) return y;		/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
+	    float u;
 	    SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
-	    y = x*x;
-	    if(y==x) return y; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u*u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hx>hy) {				/* x > y, x -= ulp */
@@ -61,15 +64,12 @@ static char rcsid[] = "$NetBSD: s_nextaf
 	if(hy>=0x7f800000) {
 	  x = x+x;	/* overflow  */
 	  if (FLT_EVAL_METHOD != 0)
-	    asm ("" : "=m"(x) : "m"(x));
+	    asm ("" : "+m"(x));
 	  return x;	/* overflow  */
 	}
-	if(hy<0x00800000) {		/* underflow */
-	    y = x*x;
-	    if(y!=x) {		/* raise underflow flag */
-	        SET_FLOAT_WORD(y,hx);
-		return y;
-	    }
+	if(hy<0x00800000) {
+	    float u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-128/s_nextafterl.c.jj	1999-07-14 02:09:34.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-128/s_nextafterl.c	2007-03-27 13:51:15.000000000 +0200
@@ -25,7 +25,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 
 #ifdef __STDC__
 	long double __nextafterl(long double x, long double y)
@@ -47,9 +47,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if(x==y) return y;		/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
+	    long double u;
 	    SET_LDOUBLE_WORDS64(x,hy&0x8000000000000000ULL,1);/* return +-minsubnormal */
-	    y = x*x;
-	    if(y==x) return y; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {			/* x > 0 */
 	    if(hx>hy||((hx==hy)&&(lx>ly))) {	/* x > y, x -= ulp */
@@ -70,12 +73,9 @@ static char rcsid[] = "$NetBSD: $";
 	}
 	hy = hx&0x7fff000000000000LL;
 	if(hy==0x7fff000000000000LL) return x+x;/* overflow  */
-	if(hy==0) {				/* underflow */
-	    y = x*x;
-	    if(y!=x) {		/* raise underflow flag */
-	        SET_LDOUBLE_WORDS64(y,hx,lx);
-		return y;
-	    }
+	if(hy==0) {
+	    long double u = x*x;		/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_LDOUBLE_WORDS64(x,hx,lx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-128/s_nexttoward.c.jj	2003-12-07 22:21:29.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-128/s_nexttoward.c	2007-03-27 14:06:31.000000000 +0200
@@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 #include <float.h>
 
 #ifdef __STDC__
@@ -53,10 +53,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
-	    double x2;
+	    double u;
 	    INSERT_WORDS(x,(u_int32_t)((hy>>32)&0x80000000),1);/* return +-minsub */
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if (hy<0||(ix>>20)>(iy>>48)-0x3c00
@@ -87,16 +89,13 @@ static char rcsid[] = "$NetBSD: $";
 	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));
+	    /* Force conversion to double.  */
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00100000) {		/* underflow */
-	    double x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        INSERT_WORDS(x2,hx,lx);
-		return x2;
-	    }
+	if(hy<0x00100000) {
+	    double u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	INSERT_WORDS(x,hx,lx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c.jj	2006-01-21 20:43:12.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-opt/s_nexttowardfd.c	2007-03-27 14:12:02.000000000 +0200
@@ -20,7 +20,10 @@
  *   Special cases:
  */
 
+#include <math.h>
+#include <math_private.h>
 #include <math_ldbl_opt.h>
+#include <float.h>
 
 float __nldbl_nexttowardf(float x, double y);
 
@@ -39,10 +42,12 @@ float __nldbl_nexttowardf(float x, doubl
 	   return x+y;
 	if((double) x==y) return y;		/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
-	    float x2;
+	    float u;
 	    SET_FLOAT_WORD(x,(u_int32_t)(hy&0x80000000)|1);/* return +-minsub*/
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x; /* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hy<0||(ix>>23)>(iy>>20)-0x380
@@ -60,13 +65,16 @@ float __nldbl_nexttowardf(float x, doubl
 		hx += 1;
 	}
 	hy = hx&0x7f800000;
-	if(hy>=0x7f800000) return x+x;	/* overflow  */
-	if(hy<0x00800000) {		/* underflow */
-	    float x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-		SET_FLOAT_WORD(x2,hx);
-		return x2;
-	    }
+	if(hy>=0x7f800000) {
+	  x = x+x;	/* overflow  */
+	  if (FLT_EVAL_METHOD != 0)
+	    /* Force conversion to float.  */
+	    asm ("" : "+m"(x));
+	  return x;
+	}
+	if(hy<0x00800000) {
+	    float u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-96/s_nextafterl.c.jj	2001-10-16 13:02:24.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-96/s_nextafterl.c	2007-03-27 13:51:30.000000000 +0200
@@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 
 #ifdef __STDC__
 	long double __nextafterl(long double x, long double y)
@@ -48,9 +48,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if(x==y) return y;		/* x=y, return y */
 	if((ix|hx|lx)==0) {			/* x == 0 */
+	    long double u;
 	    SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */
-	    y = x*x;
-	    if(y==x) return y; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(esx<0x8000) {			/* x > 0 */
 	    if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
@@ -85,13 +88,10 @@ static char rcsid[] = "$NetBSD: $";
 	    }
 	}
 	esy = esx&0x7fff;
-	if(esy==0x7fff) return x+x;	/* overflow  */
-	if(esy==0) {			/* underflow */
-	    y = x*x;
-	    if(y!=x) {		/* raise underflow flag */
-	        SET_LDOUBLE_WORDS(y,esx,hx,lx);
-		return y;
-	    }
+	if(esy==0x7fff) return x+x;		/* overflow  */
+	if(esy==0) {
+	    long double u = x*x;		/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_LDOUBLE_WORDS(x,esx,hx,lx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-96/s_nexttoward.c.jj	2006-01-14 13:08:51.000000000 +0100
+++ libc/sysdeps/ieee754/ldbl-96/s_nexttoward.c	2007-03-27 14:06:48.000000000 +0200
@@ -26,7 +26,7 @@ static char rcsid[] = "$NetBSD: $";
  */
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
 #include <float.h>
 
 #ifdef __STDC__
@@ -50,10 +50,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if((ix|lx)==0) {			/* x == 0 */
-	    double x2;
+	    double u;
 	    INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00
@@ -84,16 +86,13 @@ static char rcsid[] = "$NetBSD: $";
 	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));
+	    /* Force conversion to double.  */
+	    asm ("" : "+m"(x));
 	  return x;
 	}
-	if(hy<0x00100000) {		/* underflow */
-	    double x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        INSERT_WORDS(x2,hx,lx);
-		return x2;
-	    }
+	if(hy<0x00100000) {
+	    double u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	INSERT_WORDS(x,hx,lx);
 	return x;
--- libc/sysdeps/ieee754/ldbl-96/s_nexttowardf.c.jj	2001-06-16 05:34:42.000000000 +0200
+++ libc/sysdeps/ieee754/ldbl-96/s_nexttowardf.c	2007-03-27 14:10:56.000000000 +0200
@@ -18,7 +18,8 @@ static char rcsid[] = "$NetBSD: $";
 #endif
 
 #include "math.h"
-#include "math_private.h"
+#include <math_private.h>
+#include <float.h>
 
 #ifdef __STDC__
 	float __nexttowardf(float x, long double y)
@@ -41,10 +42,12 @@ static char rcsid[] = "$NetBSD: $";
 	   return x+y;
 	if((long double) x==y) return y;	/* x=y, return y */
 	if(ix==0) {				/* x == 0 */
-	    float x2;
+	    float u;
 	    SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
-	    x2 = x*x;
-	    if(x2==x) return x2; else return x;	/* raise underflow flag */
+	    u = math_opt_barrier (x);
+	    u = u * u;
+	    math_force_eval (u);		/* raise underflow flag */
+	    return x;
 	}
 	if(hx>=0) {				/* x > 0 */
 	    if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80
@@ -64,13 +67,16 @@ static char rcsid[] = "$NetBSD: $";
 	    }
 	}
 	hy = hx&0x7f800000;
-	if(hy>=0x7f800000) return x+x;	/* overflow  */
-	if(hy<0x00800000) {		/* underflow */
-	    float x2 = x*x;
-	    if(x2!=x) {		/* raise underflow flag */
-	        SET_FLOAT_WORD(x2,hx);
-		return x2;
-	    }
+	if(hy>=0x7f800000) {
+	  x = x+x;	/* overflow  */
+	  if (FLT_EVAL_METHOD != 0)
+	    /* Force conversion to float.  */
+	    asm ("" : "+m"(x));
+	  return x;
+	}
+	if(hy<0x00800000) {
+	    float u = x*x;			/* underflow */
+	    math_force_eval (u);		/* raise underflow flag */
 	}
 	SET_FLOAT_WORD(x,hx);
 	return x;
--- libc/sysdeps/x86_64/fpu/math_private.h.jj	2007-03-23 14:30:22.000000000 +0100
+++ libc/sysdeps/x86_64/fpu/math_private.h	2007-03-27 14:01:30.000000000 +0200
@@ -0,0 +1,21 @@
+#ifndef _MATH_PRIVATE_H
+
+#define math_opt_barrier(x) \
+({ __typeof(x) __x;					\
+   if (sizeof (x) <= sizeof (double))			\
+     __asm ("" : "=x" (__x) : "0" (x));			\
+   else							\
+     __asm ("" : "=t" (__x) : "0" (x));			\
+   __x; })
+#define math_force_eval(x) \
+do							\
+  {							\
+    if (sizeof (x) <= sizeof (double))			\
+      __asm __volatile ("" : : "x" (x));		\
+    else						\
+      __asm __volatile ("" : : "f" (x));		\
+  }							\
+while (0)
+
+#include <math/math_private.h>
+#endif

	Jakub

             reply	other threads:[~2007-03-27 13:08 UTC|newest]

Thread overview: 2+ messages / expand[flat|nested]  mbox.gz  Atom feed  top
2007-03-27 13:08 Jakub Jelinek [this message]
2007-04-01 17:07 ` [PATCH] Fix next{after,toward}l on ppc* Jakub Jelinek

Reply instructions:

You may reply publicly to this message via plain-text email
using any one of the following methods:

* Save the following mbox file, import it into your mail client,
  and reply-to-all from there: mbox

  Avoid top-posting and favor interleaved quoting:
  https://en.wikipedia.org/wiki/Posting_style#Interleaved_style

* Reply using the --to, --cc, and --in-reply-to
  switches of git-send-email(1):

  git send-email \
    --in-reply-to=20070327131242.GD1826@sunsite.mff.cuni.cz \
    --to=jakub@redhat.com \
    --cc=drepper@redhat.com \
    --cc=libc-hacker@sources.redhat.com \
    /path/to/YOUR_REPLY

  https://kernel.org/pub/software/scm/git/docs/git-send-email.html

* If your mail client supports setting the In-Reply-To header
  via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).