Use __ldexp_exp() to simplify things and improve accuracy for x near

the overflow threshold.
This commit is contained in:
das 2011-10-21 06:28:47 +00:00
parent 791a3ff0bf
commit f66ae96060
4 changed files with 10 additions and 27 deletions

View File

@ -45,7 +45,6 @@ __ieee754_cosh(double x)
{ {
double t,w; double t,w;
int32_t ix; int32_t ix;
u_int32_t lx;
/* High word of |x|. */ /* High word of |x|. */
GET_HIGH_WORD(ix,x); GET_HIGH_WORD(ix,x);
@ -72,13 +71,8 @@ __ieee754_cosh(double x)
if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x)); if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */ /* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD(lx,x); if (ix<=0x408633CE)
if (ix<0x408633CE || return __ldexp_exp(fabs(x), -1);
((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
w = __ieee754_exp(half*fabs(x));
t = half*w;
return t*w;
}
/* |x| > overflowthresold, cosh(x) overflow */ /* |x| > overflowthresold, cosh(x) overflow */
return huge*huge; return huge*huge;

View File

@ -51,11 +51,8 @@ __ieee754_coshf(float x)
if (ix < 0x42b17217) return half*__ieee754_expf(fabsf(x)); if (ix < 0x42b17217) return half*__ieee754_expf(fabsf(x));
/* |x| in [log(maxfloat), overflowthresold] */ /* |x| in [log(maxfloat), overflowthresold] */
if (ix<=0x42b2d4fc) { if (ix<=0x42b2d4fc)
w = __ieee754_expf(half*fabsf(x)); return __ldexp_expf(fabsf(x), -1);
t = half*w;
return t*w;
}
/* |x| > overflowthresold, cosh(x) overflow */ /* |x| > overflowthresold, cosh(x) overflow */
return huge*huge; return huge*huge;

View File

@ -40,9 +40,8 @@ static const double one = 1.0, shuge = 1.0e307;
double double
__ieee754_sinh(double x) __ieee754_sinh(double x)
{ {
double t,w,h; double t,h;
int32_t ix,jx; int32_t ix,jx;
u_int32_t lx;
/* High word of |x|. */ /* High word of |x|. */
GET_HIGH_WORD(jx,x); GET_HIGH_WORD(jx,x);
@ -66,12 +65,8 @@ __ieee754_sinh(double x)
if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x)); if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */ /* |x| in [log(maxdouble), overflowthresold] */
GET_LOW_WORD(lx,x); if (ix<=0x408633CE)
if (ix<0x408633CE || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { return h*2.0*__ldexp_exp(fabs(x), -1);
w = __ieee754_exp(0.5*fabs(x));
t = h*w;
return t*w;
}
/* |x| > overflowthresold, sinh(x) overflow */ /* |x| > overflowthresold, sinh(x) overflow */
return x*shuge; return x*shuge;

View File

@ -24,7 +24,7 @@ static const float one = 1.0, shuge = 1.0e37;
float float
__ieee754_sinhf(float x) __ieee754_sinhf(float x)
{ {
float t,w,h; float t,h;
int32_t ix,jx; int32_t ix,jx;
GET_FLOAT_WORD(jx,x); GET_FLOAT_WORD(jx,x);
@ -48,11 +48,8 @@ __ieee754_sinhf(float x)
if (ix < 0x42b17217) return h*__ieee754_expf(fabsf(x)); if (ix < 0x42b17217) return h*__ieee754_expf(fabsf(x));
/* |x| in [logf(maxfloat), overflowthresold] */ /* |x| in [logf(maxfloat), overflowthresold] */
if (ix<=0x42b2d4fc) { if (ix<=0x42b2d4fc)
w = __ieee754_expf((float)0.5*fabsf(x)); return h*2.0F*__ldexp_expf(fabsf(x), -1);
t = h*w;
return t*w;
}
/* |x| > overflowthresold, sinh(x) overflow */ /* |x| > overflowthresold, sinh(x) overflow */
return x*shuge; return x*shuge;