diff mbox

Fix ldbl-128ibm ilogbl near powers of 2 (bug 18029)

Message ID alpine.DEB.2.10.1502261257530.31185@digraph.polyomino.org.uk
State Committed
Headers

Commit Message

Joseph Myers Feb. 26, 2015, 12:58 p.m. UTC 
  The ldbl-128ibm implementation of ilogbl produces incorrect results
when the high part of the argument is a power of 2 and the low part a
nonzero number with the opposite sign (and so the returned exponent
should be 1 less than that of the high part).  For example, ilogbl
(0x1.ffffffffffffffp1L) returns 2 but should return 1.  (This is
similar to (fixed) bug 16740 for frexpl, and bug 18030 for logbl.)
This patch adds checks for that case.

Tested for powerpc.  Committed.

2015-02-26  Joseph Myers  <joseph@codesourcery.com>

	[BZ #18029]
	* sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c (__ieee754_ilogbl):
	Adjust exponent of power of 2 down when low part has opposite
	sign.
	* math/libm-test.inc (ilogb_test_data): Add more tests.
diff mbox

Patch

diff --git a/math/libm-test.inc b/math/libm-test.inc
index b3396b6..6045764 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7219,6 +7219,11 @@  static const struct test_f_i_data ilogb_test_data[] =
     TEST_f_i (ilogb, 1024, 10, NO_INEXACT_EXCEPTION),
     TEST_f_i (ilogb, -2000, 10, NO_INEXACT_EXCEPTION),
 
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 57
+    TEST_f_i (ilogb, 0x1.ffffffffffffffp1L, 1, NO_INEXACT_EXCEPTION),
+    TEST_f_i (ilogb, -0x1.ffffffffffffffp1L, 1, NO_INEXACT_EXCEPTION),
+#endif
+
     /* ilogb (0.0) == FP_ILOGB0 plus invalid exception  */
     TEST_f_i (ilogb, 0.0, FP_ILOGB0, NO_INEXACT_EXCEPTION|INVALID_EXCEPTION|ERRNO_EDOM),
     /* ilogb (qNaN) == FP_ILOGBNAN plus invalid exception  */
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
index aeace7c..4088238 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_ilogbl.c
@@ -31,12 +31,13 @@  static char rcsid[] = "$NetBSD: $";
 
 int __ieee754_ilogbl(long double x)
 {
-	int64_t hx;
+	int64_t hx, hxs;
 	int ix;
-	double xhi;
+	double xhi, xlo;
 
-	xhi = ldbl_high (x);
+	ldbl_unpack (x, &xhi, &xlo);
 	EXTRACT_WORDS64 (hx, xhi);
+	hxs = hx;
 	hx &= 0x7fffffffffffffffLL;
 	if(hx <= 0x0010000000000000LL) {
 	    if(hx==0)
@@ -45,7 +46,21 @@  int __ieee754_ilogbl(long double x)
 		for (ix = -1022, hx<<=11; hx>0; hx<<=1) ix -=1;
 	    return ix;
 	}
-	else if (hx<0x7ff0000000000000LL) return (hx>>52)-0x3ff;
+	else if (hx < 0x7ff0000000000000LL)
+	  {
+	    int hexp = (hx >> 52) - 0x3ff;
+	    /* If the high part is a power of 2, and the low part is
+	       nonzero with the opposite sign, the low part affects
+	       the exponent.  */
+	    if ((hx & 0x000fffffffffffffLL) == 0)
+	      {
+		int64_t lx;
+		EXTRACT_WORDS64 (lx, xlo);
+		if ((hxs ^ lx) < 0 && (lx & 0x7fffffffffffffffLL) != 0)
+		  hexp--;
+	      }
+	    return hexp;
+	  }
 	else if (FP_ILOGBNAN != INT_MAX) {
 	    /* ISO C99 requires ilogbl(+-Inf) == INT_MAX.  */
 	    if (hx==0x7ff0000000000000LL)