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, 10 deletions

diff --git a/sysdeps/ieee754/ldbl-128/s_logbl.c b/sysdeps/ieee754/ldbl-128/s_logbl.c
index 0b09b28..cf6003e 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)