Correct IBM long double frexpl.

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]
	[BZ #16619]
	* sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite.
	* math/libm-test.inc (frexp_test_data): Add tests.

diff --git a/ChangeLog b/ChangeLog
index 2a8cc57..9f70d63 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2014-04-16  Alan Modra  <amodra@gmail.com>
+
+       [BZ #16740]
+       [BZ #16619]
+       * sysdeps/ieee754/ldbl-128ibm/s_frexpl.c (__frexpl): Rewrite.
+       * math/libm-test.inc (frexp_test_data): Add tests.
+
 2014-04-15  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
        * benchtests/Makefile: Depend on libraries in build directory.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 5e6789f..a4bf0b8 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7186,6 +7186,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 (file)
--- 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);