From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 11862 invoked by alias); 27 Mar 2007 13:08:23 -0000 Received: (qmail 11835 invoked by uid 22791); 27 Mar 2007 13:08:19 -0000 X-Spam-Check-By: sourceware.org Received: from sunsite.ms.mff.cuni.cz (HELO sunsite.mff.cuni.cz) (195.113.15.26) by sourceware.org (qpsmtpd/0.31) with ESMTP; Tue, 27 Mar 2007 14:08:00 +0100 Received: from sunsite.mff.cuni.cz (localhost.localdomain [127.0.0.1]) by sunsite.mff.cuni.cz (8.13.8/8.13.8) with ESMTP id l2RDCgvQ018887; Tue, 27 Mar 2007 15:12:42 +0200 Received: (from jakub@localhost) by sunsite.mff.cuni.cz (8.13.8/8.13.8/Submit) id l2RDCga0018886; Tue, 27 Mar 2007 15:12:42 +0200 Date: Tue, 27 Mar 2007 13:08:00 -0000 From: Jakub Jelinek To: Ulrich Drepper Cc: Glibc hackers Subject: [PATCH] Fix next{after,toward}{,f,l} underflow (and overflow) handling Message-ID: <20070327131242.GD1826@sunsite.mff.cuni.cz> Reply-To: Jakub Jelinek Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.4.2.2i Mailing-List: contact libc-hacker-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-hacker-owner@sourceware.org X-SW-Source: 2007-03/txt/msg00038.txt.bz2 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 [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 -#include "math_private.h" +#include #include #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 #include +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 -#include "math_private.h" +#include #include #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 #include +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 #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 +#include #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 +#include #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 +#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 #include #include @@ -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 #include +#include #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 #include #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 #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 #include #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 +#include #include +#include 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 #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 #include #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 +#include #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 +#endif Jakub