From d7826aa149d2e85128a7ecf8fadc950ab9fe7a22 Mon Sep 17 00:00:00 2001 From: Ulrich Drepper Date: Tue, 25 Oct 2011 10:52:45 -0400 Subject: [PATCH] Use math_force_eval in more places --- ChangeLog | 22 ++++++ sysdeps/ieee754/dbl-64/e_atanh.c | 7 +- sysdeps/ieee754/dbl-64/e_j0.c | 7 +- sysdeps/ieee754/dbl-64/s_ceil.c | 34 ++++----- sysdeps/ieee754/dbl-64/s_expm1.c | 78 +++++++++----------- sysdeps/ieee754/dbl-64/s_floor.c | 33 ++++----- sysdeps/ieee754/dbl-64/s_log1p.c | 52 +++++-------- sysdeps/ieee754/dbl-64/s_round.c | 43 +++++------ sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c | 16 ++-- sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c | 16 ++-- sysdeps/ieee754/dbl-64/wordsize-64/s_round.c | 24 +++--- sysdeps/ieee754/flt-32/e_atanhf.c | 7 +- sysdeps/ieee754/flt-32/e_j0f.c | 7 +- sysdeps/ieee754/flt-32/s_ceilf.c | 14 ++-- sysdeps/ieee754/flt-32/s_expm1f.c | 52 +++++-------- sysdeps/ieee754/flt-32/s_floorf.c | 16 ++-- sysdeps/ieee754/flt-32/s_log1pf.c | 42 ++++------- sysdeps/ieee754/flt-32/s_roundf.c | 24 +++--- sysdeps/ieee754/ldbl-96/e_atanhl.c | 5 +- sysdeps/ieee754/ldbl-96/e_j0l.c | 13 ++-- sysdeps/ieee754/ldbl-96/s_roundl.c | 55 +++++++------- 21 files changed, 257 insertions(+), 310 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1046feee0b..76a21df3aa 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,27 @@ 2011-10-25 Ulrich Drepper + * sysdeps/ieee754/dbl-64/e_atanh.c: Use math_force_eval instead of a + useful if() expression. + * sysdeps/ieee754/dbl-64/e_j0.c: Likewise. + * sysdeps/ieee754/dbl-64/s_ceil.c: Likewise. + * sysdeps/ieee754/dbl-64/s_expm1.c: Likewise. + * sysdeps/ieee754/dbl-64/s_floor.c: Likewise. + * sysdeps/ieee754/dbl-64/s_log1p.c: Likewise. + * sysdeps/ieee754/dbl-64/s_round.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c: Likewise. + * sysdeps/ieee754/dbl-64/wordsize-64/s_round.c: Likewise. + * sysdeps/ieee754/flt-32/e_atanhf.c: Likewise. + * sysdeps/ieee754/flt-32/e_j0f.c: Likewise. + * sysdeps/ieee754/flt-32/s_ceilf.c: Likewise. + * sysdeps/ieee754/flt-32/s_expm1f.c: Likewise. + * sysdeps/ieee754/flt-32/s_floorf.c: Likewise. + * sysdeps/ieee754/flt-32/s_log1pf.c: Likewise. + * sysdeps/ieee754/flt-32/s_roundf.c: Likewise. + * sysdeps/ieee754/ldbl-96/e_atanhl.c: Likewise. + * sysdeps/ieee754/ldbl-96/e_j0l.c: Likewise. + * sysdeps/ieee754/ldbl-96/s_roundl.c: Likewise. + * sysdeps/x86_64/fpu/math_private.h: Use VEX encoding when possible. 2011-10-25 Andreas Schwab diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c index de3bc8f144..1f83e31981 100644 --- a/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/sysdeps/ieee754/dbl-64/e_atanh.c @@ -49,8 +49,11 @@ __ieee754_atanh (double x) double t; if (xa < 0.5) { - if (__builtin_expect (xa < 0x1.0p-28, 0) && (huge + x) > 0.0) - return x; + if (__builtin_expect (xa < 0x1.0p-28, 0)) + { + math_force_eval (huge + x); + return x; + } t = xa + xa; t = 0.5 * __log1p (t + t * xa / (1.0 - xa)); diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c index 5ebf2056bf..48584a60b4 100644 --- a/sysdeps/ieee754/dbl-64/e_j0.c +++ b/sysdeps/ieee754/dbl-64/e_j0.c @@ -111,10 +111,9 @@ __ieee754_j0(double x) return z; } if(ix<0x3f200000) { /* |x| < 2**-13 */ - if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x3e400000) return one; /* |x|<2**-27 */ - else return one - 0.25*x*x; - } + math_force_eval(huge+x); /* raise inexact if x != 0 */ + if(ix<0x3e400000) return one; /* |x|<2**-27 */ + else return one - 0.25*x*x; } z = x*x; #ifdef DO_NOT_USE_THIS diff --git a/sysdeps/ieee754/dbl-64/s_ceil.c b/sysdeps/ieee754/dbl-64/s_ceil.c index 695cae5d53..de50e29bf2 100644 --- a/sysdeps/ieee754/dbl-64/s_ceil.c +++ b/sysdeps/ieee754/dbl-64/s_ceil.c @@ -32,18 +32,17 @@ __ceil(double x) EXTRACT_WORDS(i0,i1,x); j0 = ((i0>>20)&0x7ff)-0x3ff; if(j0<20) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;i1=0;} - else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} - } + if(j0<0) { /* raise inexact if x != 0 */ + math_force_eval(huge+x); + /* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=0x80000000;i1=0;} + else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} } else { i = (0x000fffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += (0x00100000)>>j0; + i0 &= (~i); i1=0; } } else if (j0>51) { if(j0==0x400) return x+x; /* inf or NaN */ @@ -51,17 +50,16 @@ __ceil(double x) } else { i = ((u_int32_t)(0xffffffff))>>(j0-20); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) { - if(j0==20) i0+=1; - else { - j = i1 + (1<<(52-j0)); - if(j0) { + if(j0==20) i0+=1; + else { + j = i1 + (1<<(52-j0)); + if(j56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else * (vii) return 2^k(1-((E+2^-k)-r)) @@ -116,11 +112,7 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; #include "math.h" #include "math_private.h" #define one Q[0] -#ifdef __STDC__ static const double -#else -static double -#endif huge = 1.0e+300, tiny = 1.0e-300, o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */ @@ -134,12 +126,8 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ -2.01099218183624371326e-07}; /* BE8AFDB7 6E09C32D */ -#ifdef __STDC__ - double __expm1(double x) -#else - double __expm1(x) - double x; -#endif +double +__expm1(double x) { double y,hi,lo,c,t,e,hxs,hfx,r1,h2,h4,R1,R2,R3; int32_t k,xsb; @@ -153,20 +141,20 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ /* filter out huge and non-finite argument */ if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */ if(hx >= 0x40862E42) { /* if |x|>=709.78... */ - if(hx>=0x7ff00000) { + if(hx>=0x7ff00000) { u_int32_t low; GET_LOW_WORD(low,x); if(((hx&0xfffff)|low)!=0) - return x+x; /* NaN */ + return x+x; /* NaN */ else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - } - if(x > o_threshold) { + } + if(x > o_threshold) { __set_errno (ERANGE); return huge*huge; /* overflow */ } } if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */ - if(x+tiny<0.0) /* raise inexact */ + math_force_eval(x+tiny); /* raise inexact */ return tiny-one; /* return -1 */ } } @@ -187,7 +175,7 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ x = hi - lo; c = (hi-x)-lo; } - else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ + else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */ t = huge+x; /* return x with inexact flags when x!=0 */ return x - (t-(huge+x)); } @@ -212,28 +200,28 @@ Q[] = {1.0, -3.33333333333331316428e-02, /* BFA11111 111110F4 */ e -= hxs; if(k== -1) return 0.5*(x-e)-0.5; if(k==1) { - if(x < -0.25) return -2.0*(e-(x+0.5)); - else return one+2.0*(x-e); + if(x < -0.25) return -2.0*(e-(x+0.5)); + else return one+2.0*(x-e); } if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - u_int32_t high; - y = one-(e-x); + u_int32_t high; + y = one-(e-x); GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ - return y-one; + return y-one; } t = one; if(k<20) { - u_int32_t high; - SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ - y = t-(e-x); + u_int32_t high; + SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */ + y = t-(e-x); GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ } else { - u_int32_t high; + u_int32_t high; SET_HIGH_WORD(t,((0x3ff-k)<<20)); /* 2^-k */ - y = x-(e+t); - y += one; + y = x-(e+t); + y += one; GET_HIGH_WORD(high,y); SET_HIGH_WORD(y,high+(k<<20)); /* add k to y's exponent */ } diff --git a/sysdeps/ieee754/dbl-64/s_floor.c b/sysdeps/ieee754/dbl-64/s_floor.c index 5b593ca316..8b2038e69d 100644 --- a/sysdeps/ieee754/dbl-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/s_floor.c @@ -41,18 +41,16 @@ static double huge = 1.0e300; j0 = ((i0>>20)&0x7ff)-0x3ff; if(j0<20) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=i1=0;} - else if(((i0&0x7fffffff)|i1)!=0) - { i0=0xbff00000;i1=0;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=i1=0;} + else if(((i0&0x7fffffff)|i1)!=0) + { i0=0xbff00000;i1=0;} } else { i = (0x000fffff)>>j0; if(((i0&i)|i1)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00100000)>>j0; - i0 &= (~i); i1=0; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x00100000)>>j0; + i0 &= (~i); i1=0; } } else if (j0>51) { if(j0==0x400) return x+x; /* inf or NaN */ @@ -60,17 +58,16 @@ static double huge = 1.0e300; } else { i = ((u_int32_t)(0xffffffff))>>(j0-20); if((i1&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) { - if(j0==20) i0+=1; - else { - j = i1+(1<<(52-j0)); - if(jzero /* raise inexact */ - &&ax<0x3c900000) /* |x| < 2**-54 */ + math_force_eval(two54+x); /* raise inexact */ + if (ax<0x3c900000) /* |x| < 2**-54 */ return x; else return x - x*x*0.5; @@ -141,22 +125,22 @@ static double zero = 0.0; if(hx<0x43400000) { u = 1.0+x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; - c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ + k = (hu>>20)-1023; + c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */ c /= u; } else { u = x; GET_HIGH_WORD(hu,u); - k = (hu>>20)-1023; + k = (hu>>20)-1023; c = 0; } hu &= 0x000fffff; if(hu<0x6a09e) { - SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ + SET_HIGH_WORD(u,hu|0x3ff00000); /* normalize u */ } else { - k += 1; + k += 1; SET_HIGH_WORD(u,hu|0x3fe00000); /* normalize u/2 */ - hu = (0x00100000-hu)>>2; + hu = (0x00100000-hu)>>2; } f = u-1.0; } @@ -168,9 +152,9 @@ static double zero = 0.0; } R = hfsq*(1.0-0.66666666666666666*f); if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + return k*ln2_hi-((R-(k*ln2_lo+c))-f); } - s = f/(2.0+f); + s = f/(2.0+f); z = s*s; #ifdef DO_NOT_USE_THIS R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); diff --git a/sysdeps/ieee754/dbl-64/s_round.c b/sysdeps/ieee754/dbl-64/s_round.c index 94fcde0e4c..f8a38163fc 100644 --- a/sysdeps/ieee754/dbl-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1997. @@ -38,13 +38,12 @@ __round (double x) { if (j0 < 0) { - if (huge + x > 0.0) - { - i0 &= 0x80000000; - if (j0 == -1) - i0 |= 0x3ff00000; - i1 = 0; - } + math_force_eval (huge + x > 0.0); + + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3ff00000; + i1 = 0; } else { @@ -52,13 +51,12 @@ __round (double x) if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - i0 += 0x00080000 >> j0; - i0 &= ~i; - i1 = 0; - } + math_force_eval (huge + x > 0.0); + + /* Raise inexact if x != 0. */ + i0 += 0x00080000 >> j0; + i0 &= ~i; + i1 = 0; } } else if (j0 > 51) @@ -76,14 +74,13 @@ __round (double x) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - u_int32_t j = i1 + (1 << (51 - j0)); - if (j < i1) - i0 += 1; - i1 = j; - } + math_force_eval (huge + x > 0.0); + + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (51 - j0)); + if (j < i1) + i0 += 1; + i1 = j; i1 &= ~i; } diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c index e0e71558f8..346dab7995 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c @@ -32,18 +32,16 @@ __ceil(double x) EXTRACT_WORDS64(i0,x); j0 = ((i0>>52)&0x7ff)-0x3ff; if(j0<=51) { - if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=INT64_C(0x8000000000000000);} - else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} - } + if(j0<0) { /* raise inexact if x != 0 */ + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=INT64_C(0x8000000000000000);} + else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} } else { i = INT64_C(0x000fffffffffffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; + i0 &= (~i); } } else { if(j0==0x400) return x+x; /* inf or NaN */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index 8b7300bb93..5df4ab52fa 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -54,18 +54,16 @@ __floor (double x) int32_t j0 = ((i0>>52)&0x7ff)-0x3ff; if(__builtin_expect(j0<52, 1)) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=0;} - else if((i0&0x7fffffffffffffffl)!=0) - { i0=0xbff0000000000000l;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=0;} + else if((i0&0x7fffffffffffffffl)!=0) + { i0=0xbff0000000000000l;} } else { uint64_t i = (0x000fffffffffffffl)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x0010000000000000l)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x0010000000000000l)>>j0; + i0 &= (~i); } INSERT_WORDS64(x,i0); } else if (j0==0x400) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index 5bd857910c..25ff859283 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997, 2009 Free Software Foundation, Inc. + Copyright (C) 1997, 2009, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1997. @@ -37,12 +37,11 @@ __round (double x) { if (j0 < 0) { - if (huge + x > 0.0) - { - i0 &= UINT64_C(0x8000000000000000); - if (j0 == -1) - i0 |= UINT64_C(0x3ff0000000000000); - } + math_force_eval (huge + x); + + i0 &= UINT64_C(0x8000000000000000); + if (j0 == -1) + i0 |= UINT64_C(0x3ff0000000000000); } else { @@ -50,12 +49,11 @@ __round (double x) if ((i0 & i) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - i0 += UINT64_C(0x0008000000000000) >> j0; - i0 &= ~i; - } + math_force_eval (huge + x); + + /* Raise inexact if x != 0. */ + i0 += UINT64_C(0x0008000000000000) >> j0; + i0 &= ~i; } } else diff --git a/sysdeps/ieee754/flt-32/e_atanhf.c b/sysdeps/ieee754/flt-32/e_atanhf.c index ddd18ab300..d98a11ed67 100644 --- a/sysdeps/ieee754/flt-32/e_atanhf.c +++ b/sysdeps/ieee754/flt-32/e_atanhf.c @@ -49,8 +49,11 @@ __ieee754_atanhf (float x) float t; if (xa < 0.5f) { - if (__builtin_expect (xa < 0x1.0p-28f, 0) && (huge + x) > 0.0f) - return x; + if (__builtin_expect (xa < 0x1.0p-28f, 0)) + { + math_force_eval (huge + x); + return x; + } t = xa + xa; t = 0.5f * __log1pf (t + t * xa / (1.0f - xa)); diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c index d2da43f929..181c2e46d4 100644 --- a/sysdeps/ieee754/flt-32/e_j0f.c +++ b/sysdeps/ieee754/flt-32/e_j0f.c @@ -66,10 +66,9 @@ __ieee754_j0f(float x) return z; } if(ix<0x39000000) { /* |x| < 2**-13 */ - if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; - } + math_force_eval(huge+x>one); /* raise inexact if x != 0 */ + if(ix<0x32000000) return one; /* |x|<2**-27 */ + else return one - (float)0.25*x*x; } z = x*x; r = z*(R02+z*(R03+z*(R04+z*R05))); diff --git a/sysdeps/ieee754/flt-32/s_ceilf.c b/sysdeps/ieee754/flt-32/s_ceilf.c index 8a83201c15..af926600b0 100644 --- a/sysdeps/ieee754/flt-32/s_ceilf.c +++ b/sysdeps/ieee754/flt-32/s_ceilf.c @@ -29,17 +29,15 @@ __ceilf(float x) j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;} - else if(i0!=0) { i0=0x3f800000;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0<0) {i0=0x80000000;} + else if(i0!=0) { i0=0x3f800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0>0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) i0 += (0x00800000)>>j0; + i0 &= (~i); } } else { if(__builtin_expect(j0==0x80, 0)) return x+x; /* inf or NaN */ diff --git a/sysdeps/ieee754/flt-32/s_expm1f.c b/sysdeps/ieee754/flt-32/s_expm1f.c index 3f4536b906..1d707120b0 100644 --- a/sysdeps/ieee754/flt-32/s_expm1f.c +++ b/sysdeps/ieee754/flt-32/s_expm1f.c @@ -13,22 +13,14 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_expm1f.c,v 1.5 1995/05/10 20:47:11 jtc Exp $"; -#endif - #include #include "math.h" #include "math_private.h" -static const volatile float huge = 1.0e+30; -static const volatile float tiny = 1.0e-30; +static const float huge = 1.0e+30; +static const float tiny = 1.0e-30; -#ifdef __STDC__ static const float -#else -static float -#endif one = 1.0, o_threshold = 8.8721679688e+01,/* 0x42b17180 */ ln2_hi = 6.9313812256e-01,/* 0x3f317180 */ @@ -41,12 +33,8 @@ Q3 = -7.9365076090e-05, /* 0xb8a670cd */ Q4 = 4.0082177293e-06, /* 0x36867e54 */ Q5 = -2.0109921195e-07; /* 0xb457edbb */ -#ifdef __STDC__ - float __expm1f(float x) -#else - float __expm1f(x) - float x; -#endif +float +__expm1f(float x) { float y,hi,lo,c,t,e,hxs,hfx,r1; int32_t k,xsb; @@ -60,17 +48,17 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ /* filter out huge and non-finite argument */ if(hx >= 0x4195b844) { /* if |x|>=27*ln2 */ if(hx >= 0x42b17218) { /* if |x|>=88.721... */ - if(hx>0x7f800000) - return x+x; /* NaN */ + if(hx>0x7f800000) + return x+x; /* NaN */ if(hx==0x7f800000) return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */ - if(x > o_threshold) { + if(x > o_threshold) { __set_errno (ERANGE); return huge*huge; /* overflow */ } } if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */ - if(x+tiny<(float)0.0) /* raise inexact */ + math_force_eval(x+tiny);/* raise inexact */ return tiny-one; /* return -1 */ } } @@ -91,7 +79,7 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ x = hi - lo; c = (hi-x)-lo; } - else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ + else if(hx < 0x33000000) { /* when |x|<2**-25, return x */ t = huge+x; /* return x with inexact flags when x!=0 */ return x - (t-(huge+x)); } @@ -109,28 +97,28 @@ Q5 = -2.0109921195e-07; /* 0xb457edbb */ e -= hxs; if(k== -1) return (float)0.5*(x-e)-(float)0.5; if(k==1) { - if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); - else return one+(float)2.0*(x-e); + if(x < (float)-0.25) return -(float)2.0*(e-(x+(float)0.5)); + else return one+(float)2.0*(x-e); } if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */ - int32_t i; - y = one-(e-x); + int32_t i; + y = one-(e-x); GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ - return y-one; + return y-one; } t = one; if(k<23) { - int32_t i; - SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ - y = t-(e-x); + int32_t i; + SET_FLOAT_WORD(t,0x3f800000 - (0x1000000>>k)); /* t=1-2^-k */ + y = t-(e-x); GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ } else { - int32_t i; + int32_t i; SET_FLOAT_WORD(t,((0x7f-k)<<23)); /* 2^-k */ - y = x-(e+t); - y += one; + y = x-(e+t); + y += one; GET_FLOAT_WORD(i,y); SET_FLOAT_WORD(y,i+(k<<23)); /* add k to y's exponent */ } diff --git a/sysdeps/ieee754/flt-32/s_floorf.c b/sysdeps/ieee754/flt-32/s_floorf.c index dd19c6bc56..de160d2114 100644 --- a/sysdeps/ieee754/flt-32/s_floorf.c +++ b/sysdeps/ieee754/flt-32/s_floorf.c @@ -36,18 +36,16 @@ __floorf(float x) j0 = ((i0>>23)&0xff)-0x7f; if(j0<23) { if(j0<0) { /* raise inexact if x != 0 */ - if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=0;} - else if((i0&0x7fffffff)!=0) - { i0=0xbf800000;} - } + math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(i0>=0) {i0=0;} + else if((i0&0x7fffffff)!=0) + { i0=0xbf800000;} } else { i = (0x007fffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - if(huge+x>(float)0.0) { /* raise inexact flag */ - if(i0<0) i0 += (0x00800000)>>j0; - i0 &= (~i); - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) i0 += (0x00800000)>>j0; + i0 &= (~i); } } else { if(__builtin_expect(j0==0x80, 0)) return x+x; /* inf or NaN */ diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c index bd3d57635c..9e555ed570 100644 --- a/sysdeps/ieee754/flt-32/s_log1pf.c +++ b/sysdeps/ieee754/flt-32/s_log1pf.c @@ -13,18 +13,10 @@ * ==================================================== */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_log1pf.c,v 1.4 1995/05/10 20:47:48 jtc Exp $"; -#endif - #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const float -#else -static float -#endif ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ two25 = 3.355443200e+07, /* 0x4c000000 */ @@ -36,18 +28,10 @@ Lp5 = 1.8183572590e-01, /* 3E3A3325 */ Lp6 = 1.5313838422e-01, /* 3E1CD04F */ Lp7 = 1.4798198640e-01; /* 3E178897 */ -#ifdef __STDC__ static const float zero = 0.0; -#else -static float zero = 0.0; -#endif -#ifdef __STDC__ - float __log1pf(float x) -#else - float __log1pf(x) - float x; -#endif +float +__log1pf(float x) { float hfsq,f,c,s,z,R,u; int32_t k,hx,hu,ax; @@ -62,8 +46,8 @@ static float zero = 0.0; else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ } if(ax<0x31000000) { /* |x| < 2**-29 */ - if(two25+x>zero /* raise inexact */ - &&ax<0x24800000) /* |x| < 2**-54 */ + math_force_eval(two25+x); /* raise inexact */ + if (ax<0x24800000) /* |x| < 2**-54 */ return x; else return x - x*x*(float)0.5; @@ -76,37 +60,37 @@ static float zero = 0.0; if(hx<0x5a000000) { u = (float)1.0+x; GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; + k = (hu>>23)-127; /* correction term */ - c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); + c = (k>0)? (float)1.0-(u-x):x-(u-(float)1.0); c /= u; } else { u = x; GET_FLOAT_WORD(hu,u); - k = (hu>>23)-127; + k = (hu>>23)-127; c = 0; } hu &= 0x007fffff; if(hu<0x3504f7) { - SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ + SET_FLOAT_WORD(u,hu|0x3f800000);/* normalize u */ } else { - k += 1; + k += 1; SET_FLOAT_WORD(u,hu|0x3f000000); /* normalize u/2 */ - hu = (0x00800000-hu)>>2; + hu = (0x00800000-hu)>>2; } f = u-(float)1.0; } hfsq=(float)0.5*f*f; if(hu==0) { /* |f| < 2**-20 */ if(f==zero) { - if(k==0) return zero; + if(k==0) return zero; else {c += k*ln2_lo; return k*ln2_hi+c;} } R = hfsq*((float)1.0-(float)0.66666666666666666*f); if(k==0) return f-R; else - return k*ln2_hi-((R-(k*ln2_lo+c))-f); + return k*ln2_hi-((R-(k*ln2_lo+c))-f); } - s = f/((float)2.0+f); + s = f/((float)2.0+f); z = s*s; R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7)))))); if(k==0) return f-(hfsq-s*(hfsq+R)); else diff --git a/sysdeps/ieee754/flt-32/s_roundf.c b/sysdeps/ieee754/flt-32/s_roundf.c index 0b24987792..09b38cdc2a 100644 --- a/sysdeps/ieee754/flt-32/s_roundf.c +++ b/sysdeps/ieee754/flt-32/s_roundf.c @@ -1,5 +1,5 @@ /* Round float to integer away from zero. - Copyright (C) 1997 Free Software Foundation, Inc. + Copyright (C) 1997, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1997. @@ -37,12 +37,11 @@ __roundf (float x) { if (j0 < 0) { - if (huge + x > 0.0F) - { - i0 &= 0x80000000; - if (j0 == -1) - i0 |= 0x3f800000; - } + math_force_eval (huge + x > 0.0F); + + i0 &= 0x80000000; + if (j0 == -1) + i0 |= 0x3f800000; } else { @@ -50,12 +49,11 @@ __roundf (float x) if ((i0 & i) == 0) /* X is integral. */ return x; - if (huge + x > 0.0F) - { - /* Raise inexact if x != 0. */ - i0 += 0x00400000 >> j0; - i0 &= ~i; - } + math_force_eval (huge + x > 0.0F); + + /* Raise inexact if x != 0. */ + i0 += 0x00400000 >> j0; + i0 &= ~i; } } else diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c index 5a2aebef3e..0f3c7fb596 100644 --- a/sysdeps/ieee754/ldbl-96/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c @@ -52,7 +52,10 @@ __ieee754_atanhl(long double x) return (x-x)/(x-x); if(ix==0x3fff) return x/zero; - if(ix<0x3fe3&&(huge+x)>zero) return x; /* x<2**-28 */ + if(ix<0x3fe3) { + math_force_eval(huge+x); + return x; /* x<2**-28 */ + } SET_LDOUBLE_EXP(x,ix); if(ix<0x3ffe) { /* x < 0.5 */ t = x+x; diff --git a/sysdeps/ieee754/ldbl-96/e_j0l.c b/sysdeps/ieee754/ldbl-96/e_j0l.c index ce1f0f7563..abf4f109f8 100644 --- a/sysdeps/ieee754/ldbl-96/e_j0l.c +++ b/sysdeps/ieee754/ldbl-96/e_j0l.c @@ -144,13 +144,12 @@ __ieee754_j0l (long double x) } if (__builtin_expect (ix < 0x3fef, 0)) /* |x| < 2**-16 */ { - if (huge + x > one) - { /* raise inexact if x != 0 */ - if (ix < 0x3fde) /* |x| < 2^-33 */ - return one; - else - return one - 0.25 * x * x; - } + /* raise inexact if x != 0 */ + math_force_eval (huge + x); + if (ix < 0x3fde) /* |x| < 2^-33 */ + return one; + else + return one - 0.25 * x * x; } z = x * x; r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4])))); diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c index f1399cc209..833ae0d786 100644 --- a/sysdeps/ieee754/ldbl-96/s_roundl.c +++ b/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -1,5 +1,5 @@ /* Round long double to integer away from zero. - Copyright (C) 1997, 2007 Free Software Foundation, Inc. + Copyright (C) 1997, 2007, 2011 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper , 1997. @@ -38,15 +38,13 @@ __roundl (long double x) { if (j0 < 0) { - if (huge + x > 0.0) + math_force_eval (huge + x); + se &= 0x8000; + i0 = i1 = 0; + if (j0 == -1) { - se &= 0x8000; - i0 = i1 = 0; - if (j0 == -1) - { - se |= 0x3fff; - i0 = 0x80000000; - } + se |= 0x3fff; + i0 = 0x80000000; } } else @@ -55,15 +53,14 @@ __roundl (long double x) if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - if (huge + x > 0.0) - { - /* Raise inexact if x != 0. */ - u_int32_t j = i0 + (0x40000000 >> j0); - if (j < i0) - se += 1; - i0 = (j & ~i) | 0x80000000; - i1 = 0; - } + + /* Raise inexact if x != 0. */ + math_force_eval (huge + x); + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + se += 1; + i0 = (j & ~i) | 0x80000000; + i1 = 0; } } else if (j0 > 62) @@ -81,22 +78,20 @@ __roundl (long double x) /* X is integral. */ return x; - if (huge + x > 0.0) + math_force_eval (huge + x); + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (62 - j0)); + if (j < i1) { - /* Raise inexact if x != 0. */ - u_int32_t j = i1 + (1 << (62 - j0)); - if (j < i1) + u_int32_t k = i0 + 1; + if (k < i0) { - u_int32_t k = i0 + 1; - if (k < i0) - { - se += 1; - k |= 0x80000000; - } - i0 = k; + se += 1; + k |= 0x80000000; } - i1 = j; + i0 = k; } + i1 = j; i1 &= ~i; }