Message ID | 20140325023118.GG18201@bubble.grove.modra.org |
---|---|
State | New |
Headers | show |
On Tue, 25 Mar 2014, 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. Cases where the high and low double differ greatly in magnitude are bug 16619 (such test inputs can't be represented in libm-test.inc, but you should see improvements in test results for clog10l; it was large errors from clog10l that drew my attention to problems with such cases).
On Tue, Mar 25, 2014 at 12:16:51PM +0000, Joseph S. Myers wrote: > On Tue, 25 Mar 2014, 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. > > Cases where the high and low double differ greatly in magnitude are bug > 16619 (such test inputs can't be represented in libm-test.inc, but you > should see improvements in test results for clog10l; it was large errors > from clog10l that drew my attention to problems with such cases). Yes, I forgot to mention the test improvement with the patch submission. 16 clog10l tests now no longer fail, and 8 tests now have only 1 ulp error.
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); >
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);