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.

2 files changed, 30 insertions, 20 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