diff options
Diffstat (limited to 'newlib/libm/common/exp2.c')
| -rw-r--r-- | newlib/libm/common/exp2.c | 18 |
1 files changed, 13 insertions, 5 deletions
diff --git a/newlib/libm/common/exp2.c b/newlib/libm/common/exp2.c index 6579e16fd..d4883d8f6 100644 --- a/newlib/libm/common/exp2.c +++ b/newlib/libm/common/exp2.c @@ -43,6 +43,13 @@ #define C5 __exp_data.exp2_poly[4] #define C6 __exp_data.exp2_poly[5] +/* Handle cases that may overflow or underflow when computing the result that + is scale*(1+TMP) without intermediate rounding. The bit representation of + scale is in SBITS, however it has a computed exponent that may have + overflown into the sign bit so that needs to be adjusted before using it as + a double. (int32_t)KI is the k used in the argument reduction and exponent + adjustment of scale, positive k here means the result may overflow and + negative k means the result may underflow. */ static inline double specialcase (double_t tmp, uint64_t sbits, uint64_t ki) { @@ -81,6 +88,7 @@ specialcase (double_t tmp, uint64_t sbits, uint64_t ki) return check_uflow (y); } +/* Top 12 bits of a double (sign and exponent bits). */ static inline uint32_t top12 (double x) { @@ -125,22 +133,22 @@ exp2 (double x) kd -= Shift; /* k/N for int k. */ r = x - kd; /* 2^(k/N) ~= scale * (1 + tail). */ - idx = 2*(ki % N); + idx = 2 * (ki % N); top = ki << (52 - EXP_TABLE_BITS); tail = asdouble (T[idx]); /* This is only a valid scale when -1023*N < k < 1024*N. */ sbits = T[idx + 1] + top; /* exp2(x) = 2^(k/N) * 2^r ~= scale + scale * (tail + 2^r - 1). */ /* Evaluation is optimized assuming superscalar pipelined execution. */ - r2 = r*r; + r2 = r * r; /* Without fma the worst case error is 0.5/N ulp larger. */ /* Worst case error is less than 0.5+0.86/N+(abs poly error * 2^53) ulp. */ #if EXP2_POLY_ORDER == 4 - tmp = tail + r*C1 + r2*C2 + r*r2*(C3 + r*C4); + tmp = tail + r * C1 + r2 * C2 + r * r2 * (C3 + r * C4); #elif EXP2_POLY_ORDER == 5 - tmp = tail + r*C1 + r2*(C2 + r*C3) + r2*r2*(C4 + r*C5); + tmp = tail + r * C1 + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5); #elif EXP2_POLY_ORDER == 6 - tmp = tail + r*C1 + r2*(0.5 + r*C3) + r2*r2*(C4 + r*C5 + r2*C6); + tmp = tail + r * C1 + r2 * (0.5 + r * C3) + r2 * r2 * (C4 + r * C5 + r2 * C6); #endif if (unlikely (abstop == 0)) return specialcase (tmp, sbits, ki); |