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

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