diff options
Diffstat (limited to 'src/math/exp2l.c')
| -rw-r--r-- | src/math/exp2l.c | 44 | 
1 files changed, 20 insertions, 24 deletions
| diff --git a/src/math/exp2l.c b/src/math/exp2l.c index 145c20fe..8fc4037c 100644 --- a/src/math/exp2l.c +++ b/src/math/exp2l.c @@ -33,13 +33,9 @@ long double exp2l(long double x)  	return exp2(x);  }  #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -  #define TBLBITS 7  #define TBLSIZE (1 << TBLBITS) -#define BIAS    (LDBL_MAX_EXP - 1) -#define EXPMASK (BIAS + LDBL_MAX_EXP) -  static const double  redux = 0x1.8p63 / TBLSIZE,  P1    = 0x1.62e42fefa39efp-1, @@ -203,29 +199,29 @@ static const double tbl[TBLSIZE * 2] = {   */  long double exp2l(long double x)  { -	union IEEEl2bits u, v; +	union ldshape u = {x}; +	int e = u.i.se & 0x7fff;  	long double r, z; -	uint32_t hx, ix, i0; +	uint32_t i0;  	union {uint32_t u; int32_t i;} k;  	/* Filter out exceptional cases. */ -	u.e = x; -	hx = u.xbits.expsign; -	ix = hx & EXPMASK; -	if (ix >= BIAS + 14) {  /* |x| >= 16384 or x is NaN */ -		if (ix == EXPMASK) { -			if (u.xbits.man == 1ULL << 63 && hx == 0xffff) /* -inf */ +	if (e >= 0x3fff + 13) {  /* |x| >= 8192 or x is NaN */ +		if (u.i.se >= 0x3fff + 14 && u.i.se >> 15 == 0) +			/* overflow */ +			return x * 0x1p16383L; +		if (e == 0x7fff)  /* -inf or -nan */ +			return -1/x; +		if (x < -16382) { +			if (x <= -16446 || x - 0x1p63 + 0x1p63 != x) +				/* underflow */ +				FORCE_EVAL((float)(-0x1p-149/x)); +			if (x <= -16446)  				return 0; -			return x; -		} -		if (x >= 16384) { -			x *= 0x1p16383L; -			return x;  		} -		if (x <= -16446) -			return 0x1p-10000L*0x1p-10000L; -	} else if (ix < BIAS - 64)  /* |x| < 0x1p-64 */ +	} else if (e < 0x3fff - 64) {  		return 1 + x; +	}  	/*  	 * Reduce x, computing z, i0, and k. The low bits of x + redux @@ -238,13 +234,13 @@ long double exp2l(long double x)  	 * We split this into k = 0xabc and i0 = 0x12 (adjusted to  	 * index into the table), then we compute z = 0x0.003456p0.  	 */ -	u.e = x + redux; -	i0 = u.bits.manl + TBLSIZE / 2; +	u.f = x + redux; +	i0 = u.i.m + TBLSIZE / 2;  	k.u = i0 / TBLSIZE * TBLSIZE;  	k.i /= TBLSIZE;  	i0 %= TBLSIZE; -	u.e -= redux; -	z = x - u.e; +	u.f -= redux; +	z = x - u.f;  	/* Compute r = exp2l(y) = exp2lt[i0] * p(z). */  	long double t_hi = tbl[2*i0]; | 
