diff options
| -rw-r--r-- | ChangeLog | 22 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_atanh.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_j0.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_ceil.c | 34 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_expm1.c | 78 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_floor.c | 33 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_log1p.c | 52 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/s_round.c | 43 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_round.c | 24 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_atanhf.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_j0f.c | 7 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_ceilf.c | 14 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_expm1f.c | 52 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_floorf.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_log1pf.c | 42 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/s_roundf.c | 24 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_atanhl.c | 5 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_j0l.c | 13 | ||||
| -rw-r--r-- | 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 <drepper@gmail.com> + * 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 <schwab@redhat.com> 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(j<i1) i0+=1; /* got a carry */ - i1 = j; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0>0) { + if(j0==20) i0+=1; + else { + j = i1 + (1<<(52-j0)); + if(j<i1) i0+=1; /* got a carry */ + i1 = j; } - i1 &= (~i); } + i1 &= (~i); } INSERT_WORDS(x,i0,i1); return x; diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c index 324354336e..589128c08c 100644 --- a/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/sysdeps/ieee754/dbl-64/s_expm1.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; -#endif - /* expm1(x) * Returns exp(x)-1, the exponential of x minus 1. * @@ -40,38 +36,38 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... * We use a special Reme algorithm on [0,0.347] to generate - * a polynomial of degree 5 in r*r to approximate R1. The + * a polynomial of degree 5 in r*r to approximate R1. The * maximum error of this polynomial approximation is bounded * by 2**-61. In other words, * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 - * where Q1 = -1.6666666666666567384E-2, - * Q2 = 3.9682539681370365873E-4, - * Q3 = -9.9206344733435987357E-6, - * Q4 = 2.5051361420808517002E-7, - * Q5 = -6.2843505682382617102E-9; - * (where z=r*r, and the values of Q1 to Q5 are listed below) + * where Q1 = -1.6666666666666567384E-2, + * Q2 = 3.9682539681370365873E-4, + * Q3 = -9.9206344733435987357E-6, + * Q4 = 2.5051361420808517002E-7, + * Q5 = -6.2843505682382617102E-9; + * (where z=r*r, and the values of Q1 to Q5 are listed below) * with error bounded by * | 5 | -61 * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 * | | * * expm1(r) = exp(r)-1 is then computed by the following - * specific way which minimize the accumulation rounding error: + * specific way which minimize the accumulation rounding error: * 2 3 * r r [ 3 - (R1 + R1*r/2) ] * expm1(r) = r + --- + --- * [--------------------] - * 2 2 [ 6 - r*(3 - R1*r/2) ] + * 2 2 [ 6 - r*(3 - R1*r/2) ] * * To compensate the error in the argument reduction, we use * expm1(r+c) = expm1(r) + c + expm1(r)*c * ~ expm1(r) + c + r*c * Thus c+r*c will be added in as the correction terms for * expm1(r+c). Now rearrange the term to avoid optimization - * screw up: - * ( 2 2 ) - * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) + * screw up: + * ( 2 2 ) + * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) - * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) + * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) * ( ) * * = r - E @@ -87,7 +83,7 @@ static char rcsid[] = "$NetBSD: s_expm1.c,v 1.8 1995/05/10 20:47:09 jtc Exp $"; * (ii) if k=0, return r-E * (iii) if k=-1, return 0.5*(r-E)-0.5 * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) - * else return 1.0+2.0*(r-E); + * else return 1.0+2.0*(r-E); * (v) if (k<-2||k>56) 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(j<i1) i0 +=1 ; /* got a carry */ - i1=j; - } + math_force_eval(huge+x); /* raise inexact flag */ + if(i0<0) { + if(j0==20) i0+=1; + else { + j = i1+(1<<(52-j0)); + if(j<i1) i0 +=1 ; /* got a carry */ + i1=j; } - i1 &= (~i); } + i1 &= (~i); } INSERT_WORDS(x,i0,i1); return x; diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index 0a9801a931..dc79a02bb3 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -13,10 +13,6 @@ for performance improvement on pipelined processors. */ -#if defined(LIBM_SCCS) && !defined(lint) -static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; -#endif - /* double log1p(double x) * * Method : @@ -34,14 +30,14 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * 2. Approximation of log1p(f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * = 2s + s*R + * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error + * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, - * 2 4 6 8 10 12 14 + * 2 4 6 8 10 12 14 * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s - * (the values of Lp1 to Lp7 are listed in the program) + * (the values of Lp1 to Lp7 are listed in the program) * and * | 2 14 | -58.45 * | Lp1*s +...+Lp7*s - R(z) | <= 2 @@ -52,7 +48,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * log1p(f) = f - (hfsq - s*(hfsq+R)). * * 3. Finally, log1p(x) = k*ln2 + log1p(f). - * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) + * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. @@ -73,7 +69,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; * to produce the hexadecimal values shown. * * Note: Assuming log() return accurate answer, the following - * algorithm can be used to compute log1p(x) to within a few ULP: + * algorithm can be used to compute log1p(x) to within a few ULP: * * u = 1+x; * if(u==1.0) return x ; else @@ -85,11 +81,7 @@ static char rcsid[] = "$NetBSD: s_log1p.c,v 1.8 1995/05/10 20:47:46 jtc Exp $"; #include "math.h" #include "math_private.h" -#ifdef __STDC__ static const double -#else -static double -#endif ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ @@ -101,18 +93,10 @@ Lp[] = {0.0, 6.666666666666735130e-01, /* 3FE55555 55555593 */ 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 1.479819860511658591e-01}; /* 3FC2F112 DF3E5244 */ -#ifdef __STDC__ static const double zero = 0.0; -#else -static double zero = 0.0; -#endif -#ifdef __STDC__ - double __log1p(double x) -#else - double __log1p(x) - double x; -#endif +double +__log1p(double x) { double hfsq,f,c,s,z,R,u,z2,z4,z6,R1,R2,R3,R4; int32_t k,hx,hu,ax; @@ -127,8 +111,8 @@ static double zero = 0.0; else return (x-x)/(x-x); /* log1p(x<-1)=NaN */ } if(ax<0x3e200000) { /* |x| < 2**-29 */ - if(two54+x>zero /* 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 <drepper@cygnus.com>, 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 <drepper@cygnus.com>, 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 <errno.h> #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 <drepper@cygnus.com>, 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 <drepper@cygnus.com>, 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; } |