commit aa5f0ff11ad2cc85277c64cf65c723a9664e1149

Author: Alan Modra <amodra@gmail.com>

Date:   Wed Apr 16 19:33:32 2014 +0930

    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 glibc-2.17-c758a686/math/libm-test.inc glibc-2.17-c758a686/math/libm-test.inc

index 5e6789f..a4bf0b8 100644

--- glibc-2.17-c758a686/math/libm-test.inc

+++ glibc-2.17-c758a686/math/libm-test.inc

@@ -5691,6 +5691,15 @@ frexp_test (void)

   TEST_fI_f1 (frexp, 12.8L, 0.8L, 4);

   TEST_fI_f1 (frexp, -27.34L, -0.854375L, 5);

+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 106

+  TEST_fI_f1 (frexp, 1.0L-0x1p-106L, 1.0L-0x1p-106L, 0),

+  TEST_fI_f1 (frexp, 1.0L, 0.5L, 1),

+  TEST_fI_f1 (frexp, 1.0L+0x1p-105L, 0.5L+0x1p-106L, 1),

+  TEST_fI_f1 (frexp, -1.0L+0x1p-106L, -1.0L+0x1p-106L, 0),

+  TEST_fI_f1 (frexp, -1.0L, -0.5L, 1),

+  TEST_fI_f1 (frexp, -1.0L-0x1p-105L, -0.5L-0x1p-106L, 1),

+#endif

+

   END (frexp);

 }

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c	2015-06-16 15:37:44.293960271 -0400

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/s_frexpl.c	2015-06-16 13:53:25.578700428 -0400

@@ -31,57 +31,115 @@

 #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;

-

-	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;

+  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;

+

+  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 += 1LL << 52;

+	  expon -= 1;

+	}

+

+      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;

+	  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 -= 1LL << 52;

+		  expon += 1;

+		}

+	    }

+	  explo = 0;

 	}

-	*eptr += (ix>>52)-1022;

+      lx |= (explo << 52) | ixl;

+    }

-	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;

-

-	    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;

-

-	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);