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.

6 files changed, 111 insertions, 72 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c
index 2382fbb414..baa35e14d8 100644
--- a/sysdeps/ieee754/dbl-64/s_logb.c
+++ b/sysdeps/ieee754/dbl-64/s_logb.c
@@ -10,10 +10,6 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $";
-#endif
-
 /*
  * double logb(x)
  * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
@@ -23,20 +19,29 @@ static char rcsid[] = "$NetBSD: s_logb.c,v 1.8 1995/05/10 20:47:50 jtc Exp $";
 #include <math.h>
 #include <math_private.h>
 
-double __logb(double x)
+double
+__logb (double x)
 {
-	int32_t lx,ix;
-	EXTRACT_WORDS(ix,lx,x);
-	ix &= 0x7fffffff;			/* high |x| */
-	if((ix|lx)==0) return -1.0/fabs(x);
-	if(ix>=0x7ff00000) return x*x;
-	if((ix>>=20)==0) 			/* IEEE 754 logb */
-		return -1022.0;
-	else
-		return (double) (ix-1023);
+  int32_t lx, ix, rix;
+
+  EXTRACT_WORDS (ix, lx, x);
+  ix &= 0x7fffffff;		/* high |x| */
+  if ((ix | lx) == 0)
+    return -1.0 / fabs (x);
+  if (ix >= 0x7ff00000)
+    return x * x;
+  if (__builtin_expect ((rix = ix >> 20) == 0, 0))
+    {
+      /* POSIX specifies that denormal number is treated as
+         though it were normalized.  */
+      int m1 = (ix == 0) ? 0 : __builtin_clz (ix);
+      int m2 = (lx == 0) ? 0 : __builtin_clz (lx);
+      int ma = (m1 == 0) ? m2 + 32 : m1;
+      return -1022.0 + (double)(11 - ma);
+    }
+  return (double) (rix - 1023);
 }
 weak_alias (__logb, logb)
 #ifdef NO_LONG_DOUBLE
-strong_alias (__logb, __logbl)
-weak_alias (__logb, logbl)
+strong_alias (__logb, __logbl) weak_alias (__logb, logbl)
 #endif
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
index 2ad6c7ddbd..474eeef36b 100644
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
+++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c
@@ -1,5 +1,5 @@
 /* Compute radix independent exponent.
-   Copyright (C) 2011 Free Software Foundation, Inc.
+   Copyright (C) 2011, 2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
 
@@ -25,16 +25,21 @@
 double
 __logb (double x)
 {
-  int64_t ix;
+  int64_t ix, ex;
 
   EXTRACT_WORDS64 (ix, x);
   ix &= UINT64_C(0x7fffffffffffffff);
   if (ix == 0)
     return -1.0 / fabs (x);
-  unsigned int ex = ix >> 52;
+  ex = ix >> 52;
   if (ex == 0x7ff)
     return x * x;
-  return ex == 0 ? -1022.0 : (double) (ex - 1023);
+  if (__builtin_expect (ex == 0, 0))
+    {
+      int m = (ix == 0) ? 0 : __builtin_clzl (ix);
+      return -1022.0 + (double)(11 -m);
+    }
+  return (double) (ex - 1023);
 }
 weak_alias (__logb, logb)
 #ifdef NO_LONG_DOUBLE
diff --git a/sysdeps/ieee754/flt-32/s_logbf.c b/sysdeps/ieee754/flt-32/s_logbf.c
index b6aa0f057d..025c70de7e 100644
--- a/sysdeps/ieee754/flt-32/s_logbf.c
+++ b/sysdeps/ieee754/flt-32/s_logbf.c
@@ -13,23 +13,27 @@
  * ====================================================
  */
 
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_logbf.c,v 1.4 1995/05/10 20:47:51 jtc Exp $";
-#endif
-
 #include <math.h>
 #include <math_private.h>
 
-float __logbf(float x)
+float
+__logbf (float x)
 {
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;			/* high |x| */
-	if(ix==0) return (float)-1.0/fabsf(x);
-	if(ix>=0x7f800000) return x*x;
-	if((ix>>=23)==0) 			/* IEEE 754 logb */
-		return -126.0; 
-	else
-		return (float) (ix-127); 
+  int32_t ix, rix;
+
+  GET_FLOAT_WORD (ix, x);
+  ix &= 0x7fffffff;		/* high |x| */
+  if (ix == 0)
+    return (float) -1.0 / fabsf (x);
+  if (ix >= 0x7f800000)
+    return x * x;
+  if (__builtin_expect ((rix = ix >> 23) == 0, 0))
+    {
+      /* POSIX specifies that denormal number is treated as
+         though it were normalized.  */
+      int m = (ix == 0) ? 0 : __builtin_clz (ix);
+      return -126.0 + (float)(8 - m);
+    }
+  return (float) (rix - 127);
 }
 weak_alias (__logbf, logbf)
diff --git a/sysdeps/ieee754/ldbl-128/s_logbl.c b/sysdeps/ieee754/ldbl-128/s_logbl.c
index 0b09b289c2..cf6003e055 100644
--- a/sysdeps/ieee754/ldbl-128/s_logbl.c
+++ b/sysdeps/ieee754/ldbl-128/s_logbl.c
@@ -26,16 +26,27 @@ static char rcsid[] = "$NetBSD: $";
 #include <math.h>
 #include <math_private.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)==0) return -1.0/fabs(x);
-	if(hx>=0x7fff000000000000LL) return x*x;
-	if((hx>>=48)==0) 			/* IEEE 754 logb */
-		return -16382.0;
-	else
-		return (long double) (hx-0x3fff);
+  int64_t lx, hx, ex;
+
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  hx &= 0x7fffffffffffffffLL;	/* high |x| */
+  if ((hx | lx) == 0)
+    return -1.0 / fabs (x);
+  if (hx >= 0x7fff000000000000LL)
+    return x * x;
+  if ((ex = hx >> 48) == 0)	/* IEEE 754 logb */
+    {
+      /* 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 -16382.0 + (long double)(15 - ma);
+    }
+  return (long double) (ex - 16383);
 }
+
 weak_alias (__logbl, logbl)
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_logbl.c b/sysdeps/ieee754/ldbl-128ibm/s_logbl.c
index f38b129971..678b6cad57 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);
diff --git a/sysdeps/ieee754/ldbl-96/s_logbl.c b/sysdeps/ieee754/ldbl-96/s_logbl.c
index 95b644c030..d8ad4bcfcf 100644
--- a/sysdeps/ieee754/ldbl-96/s_logbl.c
+++ b/sysdeps/ieee754/ldbl-96/s_logbl.c
@@ -14,10 +14,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.h>
 #include <math_private.h>
 
-long double __logbl(long double x)
+long double
+__logbl (long double x)
 {
-	int32_t es,lx,ix;
-	GET_LDOUBLE_WORDS(es,ix,lx,x);
-	es &= 0x7fff;				/* exponent */
-	if((es|ix|lx)==0) return -1.0/fabs(x);
-	if(es==0x7fff) return x*x;
-	if(es==0) 			/* IEEE 754 logb */
-		return -16382.0;
-	else
-		return (long double) (es-0x3fff);
+  int32_t es, lx, ix;
+
+  GET_LDOUBLE_WORDS (es, ix, lx, x);
+  es &= 0x7fff;			/* exponent */
+  if ((es | ix | lx) == 0)
+    return -1.0 / fabs (x);
+  if (es == 0x7fff)
+    return x * x;
+  if (es == 0)			/* IEEE 754 logb */
+    {
+      /* POSIX specifies that denormal number is treated as
+         though it were normalized.  */
+      int m1 = (ix == 0) ? 0 : __builtin_clz (ix);
+      int m2 = (lx == 0) ? 0 : __builtin_clz (lx);
+      int ma = (m1 == 0) ? m2 + 32 : m1;
+      return -16382.0 - (long double)(ma);
+    }
+  return (long double) (es - 16383);
 }
+
 weak_alias (__logbl, logbl)