From patchwork Mon Sep 14 22:01:14 2015 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Joseph Myers X-Patchwork-Id: 517579 Return-Path: X-Original-To: incoming@patchwork.ozlabs.org Delivered-To: patchwork-incoming@bilbo.ozlabs.org Received: from sourceware.org (server1.sourceware.org [209.132.180.131]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by ozlabs.org (Postfix) with ESMTPS id EDB01140784 for ; Tue, 15 Sep 2015 08:01:32 +1000 (AEST) Authentication-Results: ozlabs.org; dkim=pass (1024-bit key; secure) header.d=sourceware.org header.i=@sourceware.org header.b=VuNQClvy; dkim-atps=neutral DomainKey-Signature: a=rsa-sha1; c=nofws; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:date:from:to:subject:message-id:mime-version :content-type; q=dns; s=default; b=a1H/sbhWl7U9XMLu8xLliuJQz2Vhq qY6HlGG6zWYpcRJDanu1+WMomGbqcNtsS2JNvzOUuuq6PyJTNNBLJ18U4N9Smd5c yiA0srXHveqJ2IOWwv8Zyq/iX5a8y5ewMQCeIO79VeZD/9stzDFkWivVFu6WuI9T aBvfyrwCH5q3OI= DKIM-Signature: v=1; a=rsa-sha1; c=relaxed; d=sourceware.org; h=list-id :list-unsubscribe:list-subscribe:list-archive:list-post :list-help:sender:date:from:to:subject:message-id:mime-version :content-type; s=default; bh=md4ct++msiKdiCBhyqY6T6z9IfY=; b=VuN QClvyqK31dKZvdIAv42KbevX+/wDpWb+OMgLzJOEeF3QZcTLUnWNmXj2ib+GTVfX lU16e1yI8oF1jWgv3rzJ9cWm84tKHvspR7i7Wg7FDqQvTrcfJ+hKCPvPJDLX0J/G 43YVMX9kkSrFx+0q0lQAefHZpVe5GcTZo1KOq+gU= Received: (qmail 126791 invoked by alias); 14 Sep 2015 22:01:25 -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 126779 invoked by uid 89); 14 Sep 2015 22:01:24 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-0.9 required=5.0 tests=AWL, BAYES_50, RCVD_IN_DNSWL_LOW, SPF_PASS autolearn=ham version=3.3.2 X-HELO: relay1.mentorg.com Date: Mon, 14 Sep 2015 22:01:14 +0000 From: Joseph Myers To: Subject: Fix exp2 missing underflows (bug 16521) [committed] Message-ID: User-Agent: Alpine 2.10 (DEB 1266 2009-07-14) MIME-Version: 1.0 Various exp2 implementations in glibc can miss underflow exceptions when the scaling down part of the calculation is exact (or, in the x86 case, when the conversion from extended precision to the target precision is exact). This patch forces the exception in a similar way to previous fixes. The x86 exp2f changes may in fact not be needed for this purpose - it's likely to be the case that no argument of type float has an exp2 result so close to an exact subnormal float value that it equals that value when rounded to 64 bits (even taking account of variation between different x86 implementations). However, they are included for consistency with the changes to exp2 and so as to fix the exp2f part of bug 18875 by ensuring that excess range and precision is removed from underflowing return values. Tested for x86_64, x86 and mips64. Committed. (auto-libm-test-out diffs omitted below.) 2015-09-14 Joseph Myers [BZ #16521] [BZ #18875] * math/e_exp2l.c (__ieee754_exp2l): Force underflow exception for small results. * sysdeps/i386/fpu/e_exp2.S (dbl_min): New object. (MO): New macro. (__ieee754_exp2): For small results, force underflow exception and remove excess range and precision from return value. * sysdeps/i386/fpu/e_exp2f.S (flt_min): New object. (MO): New macro. (__ieee754_exp2f): For small results, force underflow exception and remove excess range and precision from return value. * sysdeps/i386/fpu/e_exp2l.S (ldbl_min): New object. (MO): New macro. (__ieee754_exp2l): Force underflow exception for small results. * 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 (ldbl_min): New object. (MO): New macro. (__ieee754_exp2l): Force underflow exception for small results. * math/auto-libm-test-in: Add more tests or exp2. * math/auto-libm-test-out: Regenerated. diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in index 002a1fc..15c9146 100644 --- a/math/auto-libm-test-in +++ b/math/auto-libm-test-in @@ -1497,6 +1497,29 @@ exp2 1023 exp2 -1074 exp2 16383 exp2 -16400 +exp2 -126.125 +exp2 -126.25 +exp2 -126.375 +exp2 -126.5 +exp2 -126.625 +exp2 -126.75 +exp2 -126.875 +exp2 -1022.125 +exp2 -1022.25 +exp2 -1022.375 +exp2 -1022.5 +exp2 -1022.625 +exp2 -1022.75 +exp2 -1022.875 +exp2 -0x3.fe4e8p+8 +exp2 -0x3.fe513p+8 +exp2 -16382.125 +exp2 -16382.25 +exp2 -16382.375 +exp2 -16382.5 +exp2 -16382.625 +exp2 -16382.75 +exp2 -16382.875 exp2 0x1p-10 exp2 -0x1p-10 exp2 0x1p-20 diff --git a/math/e_exp2l.c b/math/e_exp2l.c index 8db34ad..b8cd158 100644 --- a/math/e_exp2l.c +++ b/math/e_exp2l.c @@ -38,9 +38,17 @@ __ieee754_exp2l (long double x) { int intx = (int) x; long double fractx = x - intx; + long double result; if (fabsl (fractx) < LDBL_EPSILON / 4.0L) - return __scalbnl (1.0L + fractx, intx); - return __scalbnl (__ieee754_expl (M_LN2l * fractx), intx); + result = __scalbnl (1.0L + fractx, intx); + else + result = __scalbnl (__ieee754_expl (M_LN2l * fractx), intx); + if (result < LDBL_MIN) + { + long double force_underflow = result * result; + math_force_eval (force_underflow); + } + return result; } else { diff --git a/sysdeps/i386/fpu/e_exp2.S b/sysdeps/i386/fpu/e_exp2.S index f802cf8..b75a63a 100644 --- a/sysdeps/i386/fpu/e_exp2.S +++ b/sysdeps/i386/fpu/e_exp2.S @@ -6,7 +6,24 @@ #include + .section .rodata.cst8,"aM",@progbits,8 + + .p2align 3 + .type dbl_min,@object +dbl_min: .byte 0, 0, 0, 0, 0, 0, 0x10, 0 + ASM_SIZE_DIRECTIVE(dbl_min) + +#ifdef PIC +# define MO(op) op##@GOTOFF(%ecx) +#else +# define MO(op) op +#endif + + .text ENTRY(__ieee754_exp2) +#ifdef PIC + LOAD_PIC_REG (cx) +#endif fldl 4(%esp) /* I added the following ugly construct because exp(+-Inf) resulted in NaN. The ugliness results from the bright minds at Intel. @@ -27,7 +44,22 @@ ENTRY(__ieee754_exp2) faddp /* 2^(fract(x)) */ fscale /* e^x */ fstp %st(1) - ret + fldl MO(dbl_min) + fld %st(1) + fucompp + fnstsw + sahf + jnc 3f + subl $8, %esp + cfi_adjust_cfa_offset (8) + fld %st(0) + fmul %st(0) + fstpl (%esp) + fstpl (%esp) + fldl (%esp) + addl $8, %esp + cfi_adjust_cfa_offset (-8) +3: ret 1: testl $0x200, %eax /* Test sign. */ jz 2f /* If positive, jump. */ diff --git a/sysdeps/i386/fpu/e_exp2f.S b/sysdeps/i386/fpu/e_exp2f.S index f867d0d..042c83b 100644 --- a/sysdeps/i386/fpu/e_exp2f.S +++ b/sysdeps/i386/fpu/e_exp2f.S @@ -6,7 +6,24 @@ #include + .section .rodata.cst4,"aM",@progbits,4 + + .p2align 2 + .type flt_min,@object +flt_min: .byte 0, 0, 0x80, 0 + ASM_SIZE_DIRECTIVE(flt_min) + +#ifdef PIC +# define MO(op) op##@GOTOFF(%ecx) +#else +# define MO(op) op +#endif + + .text ENTRY(__ieee754_exp2f) +#ifdef PIC + LOAD_PIC_REG (cx) +#endif flds 4(%esp) /* I added the following ugly construct because exp(+-Inf) resulted in NaN. The ugliness results from the bright minds at Intel. @@ -27,7 +44,22 @@ ENTRY(__ieee754_exp2f) faddp /* 2^(fract(x)) */ fscale /* e^x */ fstp %st(1) - ret + flds MO(flt_min) + fld %st(1) + fucompp + fnstsw + sahf + jnc 3f + subl $4, %esp + cfi_adjust_cfa_offset (4) + fld %st(0) + fmul %st(0) + fstps (%esp) + fstps (%esp) + flds (%esp) + addl $4, %esp + cfi_adjust_cfa_offset (-4) +3: ret 1: testl $0x200, %eax /* Test sign. */ jz 2f /* If positive, jump. */ diff --git a/sysdeps/i386/fpu/e_exp2l.S b/sysdeps/i386/fpu/e_exp2l.S index 2bf9a25..9a5ff95 100644 --- a/sysdeps/i386/fpu/e_exp2l.S +++ b/sysdeps/i386/fpu/e_exp2l.S @@ -6,7 +6,24 @@ #include + .section .rodata.cst16,"aM",@progbits,16 + .p2align 4 + .type ldbl_min,@object +ldbl_min: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x1, 0 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(ldbl_min) + +#ifdef PIC +# define MO(op) op##@GOTOFF(%ecx) +#else +# define MO(op) op +#endif + + .text ENTRY(__ieee754_exp2l) +#ifdef PIC + LOAD_PIC_REG (cx) +#endif fldt 4(%esp) /* I added the following ugly construct because exp(+-Inf) resulted in NaN. The ugliness results from the bright minds at Intel. @@ -35,7 +52,17 @@ ENTRY(__ieee754_exp2l) faddp /* 2^(fract(x)) */ fscale /* e^x */ fstp %st(1) - ret + /* Ensure underflow for tiny result. */ + fldt MO(ldbl_min) + fld %st(1) + fucompp + fnstsw + sahf + jnc 4f + fld %st + fmul %st + fstp %st +4: ret 1: testl $0x200, %eax /* Test sign. */ jz 2f /* If positive, jump. */ diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index 30f0a8f..01df621 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -120,7 +120,15 @@ __ieee754_exp2 (double x) if (!unsafe) return result; else - return result * scale_u.d; + { + result *= scale_u.d; + if (result < DBL_MIN) + { + double force_underflow = result * result; + math_force_eval (force_underflow); + } + return result; + } } else /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c index 0b75a7e..bc62042 100644 --- a/sysdeps/ieee754/flt-32/e_exp2f.c +++ b/sysdeps/ieee754/flt-32/e_exp2f.c @@ -109,7 +109,15 @@ __ieee754_exp2f (float x) if (!unsafe) return result; else - return result * scale_u.f; + { + result *= scale_u.f; + if (result < FLT_MIN) + { + float force_underflow = result * result; + math_force_eval (force_underflow); + } + return result; + } } /* Exceptional cases: */ else if (isless (x, himark)) diff --git a/sysdeps/x86_64/fpu/e_exp2l.S b/sysdeps/x86_64/fpu/e_exp2l.S index 7d42a93..d634ad3 100644 --- a/sysdeps/x86_64/fpu/e_exp2l.S +++ b/sysdeps/x86_64/fpu/e_exp2l.S @@ -7,6 +7,20 @@ #include + .section .rodata.cst16,"aM",@progbits,16 + .p2align 4 + .type ldbl_min,@object +ldbl_min: .byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x1, 0 + .byte 0, 0, 0, 0, 0, 0 + ASM_SIZE_DIRECTIVE(ldbl_min) + +#ifdef PIC +# define MO(op) op##(%rip) +#else +# define MO(op) op +#endif + + .text ENTRY(__ieee754_exp2l) fldt 8(%rsp) /* I added the following ugly construct because exp(+-Inf) resulted @@ -36,7 +50,16 @@ ENTRY(__ieee754_exp2l) faddp /* 2^(fract(x)) */ fscale /* e^x */ fstp %st(1) - ret + /* Ensure underflow for tiny result. */ + fldt MO(ldbl_min) + fld %st(1) + fucomip %st(1), %st + fstp %st + jnc 4f + fld %st + fmul %st + fstp %st +4: ret 1: testl $0x200, %eax /* Test sign. */ jz 2f /* If positive, jump. */