From patchwork Tue Mar 25 02:31:18 2014 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Alan Modra X-Patchwork-Id: 333217 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 61EAC140085 for ; Tue, 25 Mar 2014 13:31:39 +1100 (EST) 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=nDw5PlWgWjE8aiKDL6yg0mq0JDjGU mepV5g/YvFHX7vBmYgAMMnGQFnqDLy0k2ahJW1hF5XvMSldS4CVVZqbG9kMjqAmi qo4zbour/ovcW+FFZWbpRvQxHMTjI+nj+hQioy+1gNNRGeoBAD3SMBBkHxA9UuZO sCCN/Rg1MdbCx8= 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=QQZEgCTAq0CPVSKx67AViQad6Yk=; b=dcr x5+yiXvG1it7AK9QerCKGCVR30FudIvg9XvkBl+6ZJFSJc7dm/S0dXc3IyzM/yLC LlAhH4/kjAGKD2m9ApxvCrCxctYZaqCX3cZpKbYaavDY3LGsgAWtfpxS6lw0kFgh 7TzXsFJJY4ss1uP9j8pbf+OPrXsXDa77Q07pTvzw= Received: (qmail 9790 invoked by alias); 25 Mar 2014 02:31:33 -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 9731 invoked by uid 89); 25 Mar 2014 02:31:27 -0000 Authentication-Results: sourceware.org; auth=none X-Virus-Found: No X-Spam-SWARE-Status: No, score=-2.2 required=5.0 tests=AWL, BAYES_00, FREEMAIL_FROM, RCVD_IN_DNSWL_LOW, SPF_PASS autolearn=ham version=3.3.2 X-HELO: mail-pa0-f53.google.com X-Received: by 10.66.141.197 with SMTP id rq5mr76247695pab.64.1395714683133; Mon, 24 Mar 2014 19:31:23 -0700 (PDT) Date: Tue, 25 Mar 2014 13:01:18 +1030 From: Alan Modra To: libc-alpha@sourceware.org Subject: [PowerPC64] Correct IBM long double frexpl Message-ID: <20140325023118.GG18201@bubble.grove.modra.org> MIME-Version: 1.0 Content-Disposition: inline User-Agent: Mutt/1.5.21 (2010-09-15) Besides fixing the bugzilla, this also fixes corner-cases where the high and low double differ greatly in magnitude, and handles a denormal input without resorting to a fp rescale. [BZ #16740] * sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite. * math/libm-test.inc (frexp_test_data): Add tests. diff --git a/math/libm-test.inc b/math/libm-test.inc index e97b18a..01e09d2 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -7202,6 +7202,15 @@ static const struct test_f_f1_data frexp_test_data[] = TEST_fI_f1 (frexp, 12.8L, 0.8L, 4, NO_INEXACT_EXCEPTION), TEST_fI_f1 (frexp, -27.34L, -0.854375L, 5, NO_INEXACT_EXCEPTION), + +#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 106 + TEST_fI_f1 (frexp, 1.0L-0x1p-106L, 1.0L-0x1p-106L, 0, NO_INEXACT_EXCEPTION), + TEST_fI_f1 (frexp, 1.0L, 0.5L, 1, NO_INEXACT_EXCEPTION), + TEST_fI_f1 (frexp, 1.0L+0x1p-105L, 0.5L+0x1p-106L, 1, NO_INEXACT_EXCEPTION), + TEST_fI_f1 (frexp, -1.0L+0x1p-106L, -1.0L+0x1p-106L, 0, NO_INEXACT_EXCEPTION), + TEST_fI_f1 (frexp, -1.0L, -0.5L, 1, NO_INEXACT_EXCEPTION), + TEST_fI_f1 (frexp, -1.0L-0x1p-105L, -0.5L-0x1p-106L, 1, NO_INEXACT_EXCEPTION), +#endif }; static void diff --git a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c index 7e40663..a644f92 100644 --- a/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c @@ -31,57 +31,115 @@ static char rcsid[] = "$NetBSD: $"; #include #include -static const long double -two107 = 162259276829213363391578010288128.0; /* 0x4670000000000000, 0 */ - long double __frexpl(long double x, int *eptr) { - uint64_t hx, lx, ix, ixl; - int64_t explo; - double xhi, xlo; + uint64_t hx, lx, ix, ixl; + int64_t explo, expon; + double xhi, xlo; + + ldbl_unpack (x, &xhi, &xlo); + EXTRACT_WORDS64 (hx, xhi); + EXTRACT_WORDS64 (lx, xlo); + ixl = 0x7fffffffffffffffULL & lx; + ix = 0x7fffffffffffffffULL & hx; + expon = 0; + if (ix >= 0x7ff0000000000000ULL || ix == 0) + { + /* 0,inf,nan. */ + *eptr = expon; + return x; + } + expon = ix >> 52; + if (expon == 0) + { + /* Denormal high double, the low double must be 0.0. */ + int cnt; + + /* Normalize. */ + if (sizeof (ix) == sizeof (long)) + cnt = __builtin_clzl (ix); + else if ((ix >> 32) != 0) + cnt = __builtin_clzl ((long) (ix >> 32)); + else + cnt = __builtin_clzl ((long) ix) + 32; + cnt = cnt - 12; + expon -= cnt; + ix <<= cnt + 1; + } + expon -= 1022; + ix &= 0x000fffffffffffffULL; + hx &= 0x8000000000000000ULL; + hx |= (1022LL << 52) | ix; - ldbl_unpack (x, &xhi, &xlo); - EXTRACT_WORDS64 (hx, xhi); - EXTRACT_WORDS64 (lx, xlo); - ixl = 0x7fffffffffffffffULL&lx; - ix = 0x7fffffffffffffffULL&hx; - *eptr = 0; - if(ix>=0x7ff0000000000000ULL||ix==0) return x; /* 0,inf,nan */ - if (ix<0x0010000000000000ULL) { /* subnormal */ - x *= two107; - xhi = ldbl_high (x); - EXTRACT_WORDS64 (hx, xhi); - ix = hx&0x7fffffffffffffffULL; - *eptr = -107; + if (ixl != 0) + { + /* If the high double is an exact power of two and the low + double has the opposite sign, then the exponent calculated + from the high double is one too big. */ + if (ix == 0 + && (int64_t) (hx ^ lx) < 0) + { + hx += 1L << 52; + expon -= 1; } - *eptr += (ix>>52)-1022; - if (ixl != 0ULL) { - explo = (ixl>>52) - (ix>>52) + 0x3fe; - if ((ixl&0x7ff0000000000000ULL) == 0LL) { - /* the lower double is a denormal so we need to correct its - mantissa and perhaps its exponent. */ - int cnt; + explo = ixl >> 52; + if (explo == 0) + { + /* The low double started out as a denormal. Normalize its + mantissa and adjust the exponent. */ + int cnt; - if (sizeof (ixl) == sizeof (long)) - cnt = __builtin_clzl (ixl); - else if ((ixl >> 32) != 0) - cnt = __builtin_clzl ((long) (ixl >> 32)); - else - cnt = __builtin_clzl ((long) ixl) + 32; - cnt = cnt - 12; - lx = (lx&0x8000000000000000ULL) | ((explo-cnt)<<52) - | ((ixl<<(cnt+1))&0x000fffffffffffffULL); - } else - lx = (lx&0x800fffffffffffffULL) | (explo<<52); - } else - lx = 0ULL; + if (sizeof (ixl) == sizeof (long)) + cnt = __builtin_clzl (ixl); + else if ((ixl >> 32) != 0) + cnt = __builtin_clzl ((long) (ixl >> 32)); + else + cnt = __builtin_clzl ((long) ixl) + 32; + cnt = cnt - 12; + explo -= cnt; + ixl <<= cnt + 1; + } + + /* With variable precision we can't assume much about the + magnitude of the returned low double. It may even be a + denormal. */ + explo -= expon; + ixl &= 0x000fffffffffffffULL; + lx &= 0x8000000000000000ULL; + if (explo <= 0) + { + /* Handle denormal low double. */ + if (explo > -52) + { + ixl |= 1LL << 52; + ixl >>= 1 - explo; + } + else + { + ixl = 0; + lx = 0; + if ((hx & 0x7ff0000000000000ULL) == (1023LL << 52)) + { + /* Oops, the adjustment we made above for values a + little smaller than powers of two turned out to + be wrong since the returned low double will be + zero. This can happen if the input was + something weird like 0x1p1000 - 0x1p-1000. */ + hx -= 1L << 52; + expon += 1; + } + } + explo = 0; + } + lx |= (explo << 52) | ixl; + } - hx = (hx&0x800fffffffffffffULL) | 0x3fe0000000000000ULL; - INSERT_WORDS64 (xhi, hx); - INSERT_WORDS64 (xlo, lx); - x = ldbl_pack (xhi, xlo); - return x; + INSERT_WORDS64 (xhi, hx); + INSERT_WORDS64 (xlo, lx); + x = ldbl_pack (xhi, xlo); + *eptr = expon; + return x; } #ifdef IS_IN_libm long_double_symbol (libm, __frexpl, frexpl);