mirror/glibc
2
0
Fork 0
mirror of git://sourceware.org/git/glibc.git synced 2024-11-21 01:12:26 +08:00
Code Issues Packages Projects Releases Wiki Activity

Use math_force_eval in more places

Browse Source
This commit is contained in:
Ulrich Drepper 2011-10-25 10:52:45 -04:00
parent 31ea014d8b
commit d7826aa149
21 changed files with 257 additions and 310 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
ChangeLog
View File

@ -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>

7
sysdeps/ieee754/dbl-64/e_atanh.c
View File

@ -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));

7
sysdeps/ieee754/dbl-64/e_j0.c
View File

@ -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

34
sysdeps/ieee754/dbl-64/s_ceil.c
View File

@ -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;

78
sysdeps/ieee754/dbl-64/s_expm1.c
View File

@ -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 */
}

33
sysdeps/ieee754/dbl-64/s_floor.c
View File

@ -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;

52
sysdeps/ieee754/dbl-64/s_log1p.c
View File

@ -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))))));

43
sysdeps/ieee754/dbl-64/s_round.c
View File

@ -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;
}

16
sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c
View File

@ -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 */

16
sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c
View File

@ -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)

24
sysdeps/ieee754/dbl-64/wordsize-64/s_round.c
View File

@ -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

7
sysdeps/ieee754/flt-32/e_atanhf.c
View File

@ -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));

7
sysdeps/ieee754/flt-32/e_j0f.c
View File

@ -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)));

14
sysdeps/ieee754/flt-32/s_ceilf.c
View File

@ -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 */

52
sysdeps/ieee754/flt-32/s_expm1f.c
View File

@ -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 */
}

16
sysdeps/ieee754/flt-32/s_floorf.c
View File

@ -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 */

42
sysdeps/ieee754/flt-32/s_log1pf.c
View File

@ -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

24
sysdeps/ieee754/flt-32/s_roundf.c
View File

@ -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

5
sysdeps/ieee754/ldbl-96/e_atanhl.c
View File

@ -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;

13
sysdeps/ieee754/ldbl-96/e_j0l.c
View File

@ -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]))));

55
sysdeps/ieee754/ldbl-96/s_roundl.c
View File

@ -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;
}