Fix for logb/logbf/logbl (bugs 13954/13955/13956)

POSIX 2008 states that if the input for 'logb[f|l]' is a subnormal number
it should be treated as if it were normalized.  This means the
implementation should calculate the log2 of the mantissa and add it to the
subnormal exponent (-126 for float and -1022 for double and IBM long
double).  This patch takes care of that.

1 files changed, 21 insertions, 14 deletions

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c
index f38b129..678b6ca 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c
@@ -13,10 +13,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: $";
-#endif
-
 /*
  * long double logbl(x)
  * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
@@ -27,16 +23,27 @@ static char rcsid[] = "$NetBSD: $";
 #include <math_private.h>
 #include <math_ldbl_opt.h>
 
-long double __logbl(long double x)
+long double
+__logbl (long double x)
 {
-	int64_t lx,hx;
-	GET_LDOUBLE_WORDS64(hx,lx,x);
-	hx &= 0x7fffffffffffffffLL;		/* high |x| */
-	if((hx|(lx&0x7fffffffffffffffLL))==0) return -1.0/fabs(x);
-	if(hx>=0x7ff0000000000000LL) return x*x;
-	if((hx>>=52)==0) 			/* IEEE 754 logb */
-		return -1022.0;
-	else
-		return (long double) (hx-0x3ff);
+  int64_t lx, hx, rhx;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  hx &= 0x7fffffffffffffffLL;	/* high |x| */
+  if ((hx | (lx & 0x7fffffffffffffffLL)) == 0)
+    return -1.0 / fabs (x);
+  if (hx >= 0x7ff0000000000000LL)
+    return x * x;
+  if (__builtin_expect ((rhx = hx >> 52) == 0, 0))
+    {
+      /* POSIX specifies that denormal number is treated as
+         though it were normalized.  */
+      int m1 = (hx == 0) ? 0 : __builtin_clzll (hx);
+      int m2 = (lx == 0) ? 0 : __builtin_clzll (lx);
+      int ma = (m1 == 0) ? m2 + 64 : m1;
+      return -1022.0 + (long double)(11 - ma);
+    }
+  return (long double) (rhx - 1023);
 }
+
 long_double_symbol (libm, __logbl, logbl);