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>
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));
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
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 */
} 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;
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.
*
* = 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
* (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))
#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 */
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;
/* 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 */
}
}
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));
}
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 */
}
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 */
} 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;
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 :
* 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
* 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.
* 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
#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 */
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;
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;
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;
}
}
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))))));
/* 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.
{
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
{
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)
/* 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;
}
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 */
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)
/* 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.
{
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
{
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
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));
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)));
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 */
* ====================================================
*/
-#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 */
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;
/* 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 */
}
}
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));
}
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 */
}
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 */
* ====================================================
*/
-#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 */
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;
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;
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
/* 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.
{
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
{
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
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;
}
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]))));
/* 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.
{
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
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)
/* 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;
}
projects / platform / upstream / glibc.git / commitdiff