This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Re: [PowerPC64] Correct IBM long double frexpl
- From: Adhemerval Zanella <azanella at linux dot vnet dot ibm dot com>
- To: libc-alpha at sourceware dot org
- Date: Wed, 02 Apr 2014 09:37:35 -0300
- Subject: Re: [PowerPC64] Correct IBM long double frexpl
- Authentication-results: sourceware.org; auth=none
- References: <20140325023118 dot GG18201 at bubble dot grove dot modra dot org>
Hi Alan,
LGTM, thanks.
On 24-03-2014 23:31, Alan Modra wrote:
> 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 <math_private.h>
> #include <math_ldbl_opt.h>
>
> -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);
>