Fixed fdlibm[+cygnus] logbf() and logb() on denormals. Adjustment

according to the highest nonzero bit in a denormal was missing.

fdlibm ilogbf() and ilogb() have always had the adjustment, but only
use a small part of their method for handling denormals; use the
normalization method in log[f]() for the main part.

2 files changed, 20 insertions, 8 deletions

diff --git a/lib/msun/src/s_logb.c b/lib/msun/src/s_logb.c
index 5c6ea6c..f667bf2 100644
--- a/lib/msun/src/s_logb.c
+++ b/lib/msun/src/s_logb.c
@@ -23,6 +23,9 @@ static char rcsid[] = "$FreeBSD$";
 #include "math.h"
 #include "math_private.h"
 
+static const double
+two54 = 1.80143985094819840000e+16;	/* 43500000 00000000 */
+
 double
 logb(double x)
 {
@@ -31,8 +34,11 @@ logb(double 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);
+	if(ix<0x00100000) {
+		x *= two54;		 /* convert subnormal x to normal */
+		GET_FLOAT_WORD(ix,x);
+		ix &= 0x7fffffff;
+		return (float) ((ix>>20)-1023-54);
+	} else
+		return (double) ((ix>>20)-1023);
 }
diff --git a/lib/msun/src/s_logbf.c b/lib/msun/src/s_logbf.c
index e673d96..97ccd02 100644
--- a/lib/msun/src/s_logbf.c
+++ b/lib/msun/src/s_logbf.c
@@ -20,6 +20,9 @@ static char rcsid[] = "$FreeBSD$";
 #include "math.h"
 #include "math_private.h"
 
+static const float
+two25 = 3.355443200e+07;		/* 0x4c000000 */
+
 float
 logbf(float x)
 {
@@ -28,8 +31,11 @@ logbf(float 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);
+	if(ix<0x00800000) {
+		x *= two25;		 /* convert subnormal x to normal */
+		GET_FLOAT_WORD(ix,x);
+		ix &= 0x7fffffff;
+		return (float) ((ix>>23)-127-25);
+	} else
+		return (float) ((ix>>23)-127);
 }