From patchwork Thu Feb 12 19:03:38 2015 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Joseph Myers X-Patchwork-Id: 5057 Received: (qmail 16286 invoked by alias); 12 Feb 2015 19:03:47 -0000 Mailing-List: contact libc-alpha-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Unsubscribe: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: libc-alpha-owner@sourceware.org Delivered-To: mailing list libc-alpha@sourceware.org Received: (qmail 16266 invoked by uid 89); 12 Feb 2015 19:03:46 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-2.3 required=5.0 tests=AWL, BAYES_00, RCVD_IN_DNSWL_LOW, SPF_PASS autolearn=ham version=3.3.2 X-HELO: relay1.mentorg.com Date: Thu, 12 Feb 2015 19:03:38 +0000 From: Joseph Myers To: Subject: Fix exp2 spurious underflows (bug 16560) [committed] Message-ID: User-Agent: Alpine 2.10 (DEB 1266 2009-07-14) MIME-Version: 1.0 This patch fixes the remaining part of bug 16560, spurious underflows from exp2 of arguments close to 0 (when the result is close to 1, so should not underflow), by just using 1+x instead of a more complicated calculation when the argument is sufficiently small. Tested for x86_64, x86 and mips64. Committed. (auto-libm-test-out diffs omitted below.) 2015-02-12 Joseph Myers [BZ #16560] * math/e_exp2l.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine and redefine. (__ieee754_exp2l): Do not multiply small fractional parts by M_LN2l. * sysdeps/i386/fpu/e_exp2l.S (__ieee754_exp2l): Just add 1 to small argument. * sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Likewise. * sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise. * sysdeps/x86_64/fpu/e_exp2l.S (__ieee754_exp2l): Likewise. * math/auto-libm-test-in: Add more tests of exp2. * math/auto-libm-test-out: Regenerated. diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index ad43113..6bcfd54 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -938,6 +938,24 @@ exp2 1023 exp2 -1074 exp2 16383 exp2 -16400 +exp2 0x1p-10 +exp2 -0x1p-10 +exp2 0x1p-20 +exp2 -0x1p-20 +exp2 0x1p-30 +exp2 -0x1p-30 +exp2 0x1p-40 +exp2 -0x1p-40 +exp2 0x1p-50 +exp2 -0x1p-50 +exp2 0x1p-60 +exp2 -0x1p-60 +exp2 0x1p-100 +exp2 -0x1p-100 +exp2 min +exp2 -min +exp2 min_subnorm +exp2 -min_subnorm expm1 0 expm1 -0 diff --git a/math/e_exp2l.c b/math/e_exp2l.c index bb7feef..9eb7baf 100644 --- a/math/e_exp2l.c +++ b/math/e_exp2l.c @@ -20,6 +20,13 @@ #include #include +/* To avoid spurious underflows, use this definition to treat IBM long + double as approximating an IEEE-style format. */ +#if LDBL_MANT_DIG == 106 +# undef LDBL_EPSILON +# define LDBL_EPSILON 0x1p-106L +#endif + long double __ieee754_exp2l (long double x) { @@ -31,6 +38,8 @@ __ieee754_exp2l (long double x) { int intx = (int) x; long double fractx = x - intx; + if (fabsl (fractx) < LDBL_EPSILON / 4.0L) + return __scalbnl (1.0L + fractx, intx); return __scalbnl (__ieee754_expl (M_LN2l * fractx), intx); } else diff --git a/sysdeps/i386/fpu/e_exp2l.S b/sysdeps/i386/fpu/e_exp2l.S index 203dd00..2bf9a25 100644 --- a/sysdeps/i386/fpu/e_exp2l.S +++ b/sysdeps/i386/fpu/e_exp2l.S @@ -18,7 +18,15 @@ ENTRY(__ieee754_exp2l) andb %ah, %dh cmpb $0x05, %dh je 1f /* Is +-Inf, jump. */ - fld %st + movzwl 4+8(%esp), %eax + andl $0x7fff, %eax + cmpl $0x3fbe, %eax + jge 3f + /* Argument's exponent below -65, result rounds to 1. */ + fld1 + faddp + ret +3: fld %st frndint /* int(x) */ fsubr %st,%st(1) /* fract(x) */ fxch diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index 3666c6a..f964a5a 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -61,6 +61,9 @@ __ieee754_exp2 (double x) double rx, x22, result; union ieee754_double ex2_u, scale_u; + if (fabs (x) < DBL_EPSILON / 4.0) + return 1.0 + x; + { SET_RESTORE_ROUND_NOEX (FE_TONEAREST); diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 01cd444..f3e3a8e 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -54,6 +54,9 @@ __ieee754_exp2f (float x) float rx, x22, result; union ieee754_float ex2_u, scale_u; + if (fabsf (x) < FLT_EPSILON / 4.0f) + return 1.0f + x; + { SET_RESTORE_ROUND_NOEXF (FE_TONEAREST); diff --git a/sysdeps/x86_64/fpu/e_exp2l.S b/sysdeps/x86_64/fpu/e_exp2l.S index 7abf425..7d42a93 100644 --- a/sysdeps/x86_64/fpu/e_exp2l.S +++ b/sysdeps/x86_64/fpu/e_exp2l.S @@ -19,7 +19,15 @@ ENTRY(__ieee754_exp2l) andb %ah, %dh cmpb $0x05, %dh je 1f /* Is +-Inf, jump. */ - fld %st + movzwl 8+8(%rsp), %eax + andl $0x7fff, %eax + cmpl $0x3fbe, %eax + jge 3f + /* Argument's exponent below -65, result rounds to 1. */ + fld1 + faddp + ret +3: fld %st frndint /* int(x) */ fsubr %st,%st(1) /* fract(x) */ fxch