From aa5f0ff11ad2cc85277c64cf65c723a9664e1149 Mon Sep 17 00:00:00 2001
From: Alan Modra <amodra@gmail.com>
Date: Wed, 16 Apr 2014 19:33:32 +0930
Subject: 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.
---
 sysdeps/ieee754/ldbl-128ibm/s_frexpl.c | 146 +++++++++++++++++++++++----------
 1 file changed, 102 insertions(+), 44 deletions(-)

(limited to 'sysdeps/ieee754/ldbl-128ibm/s_frexpl.c')

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);
-- 
cgit v1.1