* libc/include/machine/ieeefp.h: Comment about new configuration

	macros _FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS.
	* libm/common/fdlib.h: Define new macros for testing floats.
	* libm/common/sf_*: Use them.
	* libm/math/ef_*: Likewise.
	* libm/math/sf_*: Likewise.

11 files changed, 158 insertions, 47 deletions

diff --git a/newlib/libm/common/fdlibm.h b/newlib/libm/common/fdlibm.h
index 752d24246..5ca1a917b 100644
--- a/newlib/libm/common/fdlibm.h
+++ b/newlib/libm/common/fdlibm.h
@@ -18,14 +18,116 @@
 /* CYGNUS LOCAL: Default to XOPEN_MODE.  */
 #define _XOPEN_MODE
 
+/* Most routines need to check whether a float is finite, infinite, or not a
+   number, and many need to know whether the result of an operation will
+   overflow.  These conditions depend on whether the largest exponent is
+   used for NaNs & infinities, or whether it's used for finite numbers.  The
+   macros below wrap up that kind of information:
+
+   FLT_UWORD_IS_FINITE(X)
+	True if a positive float with bitmask X is finite.
+
+   FLT_UWORD_IS_NAN(X)
+	True if a positive float with bitmask X is not a number.
+
+   FLT_UWORD_IS_INFINITE(X)
+	True if a positive float with bitmask X is +infinity.
+
+   FLT_UWORD_MAX
+	The bitmask of FLT_MAX.
+
+   FLT_UWORD_HALF_MAX
+	The bitmask of FLT_MAX/2.
+
+   FLT_UWORD_EXP_MAX
+	The bitmask of the largest finite exponent (129 if the largest
+	exponent is used for finite numbers, 128 otherwise).
+
+   FLT_UWORD_LOG_MAX
+	The bitmask of log(FLT_MAX), rounded down.  This value is the largest
+	input that can be passed to exp() without producing overflow.
+
+   FLT_UWORD_LOG_2MAX
+	The bitmask of log(2*FLT_MAX), rounded down.  This value is the
+	largest input than can be passed to cosh() without producing
+	overflow.
+
+   FLT_LARGEST_EXP
+	The largest biased exponent that can be used for finite numbers
+	(255 if the largest exponent is used for finite numbers, 254
+	otherwise) */
+
+#ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
+#define FLT_UWORD_IS_FINITE(x) 1
+#define FLT_UWORD_IS_NAN(x) 0
+#define FLT_UWORD_IS_INFINITE(x) 0
+#define FLT_UWORD_MAX 0x7fffffff
+#define FLT_UWORD_EXP_MAX 0x43010000
+#define FLT_UWORD_LOG_MAX 0x42b2d4fc
+#define FLT_UWORD_LOG_2MAX 0x42b437e0
+#define HUGE ((float)0X1.FFFFFEP128)
+#else
+#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
+#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
+#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
+#define FLT_UWORD_MAX 0x7f7fffff
+#define FLT_UWORD_EXP_MAX 0x43000000
+#define FLT_UWORD_LOG_MAX 0x42b17217
+#define FLT_UWORD_LOG_2MAX 0x42b2d4fc
+#define HUGE ((float)3.40282346638528860e+38)
+#endif
+#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1<<23))
+#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
+
+/* Many routines check for zero and subnormal numbers.  Such things depend
+   on whether the target supports denormals or not:
+
+   FLT_UWORD_IS_ZERO(X)
+	True if a positive float with bitmask X is +0.	Without denormals,
+	any float with a zero exponent is a +0 representation.	With
+	denormals, the only +0 representation is a 0 bitmask.
+
+   FLT_UWORD_IS_SUBNORMAL(X)
+	True if a non-zero positive float with bitmask X is subnormal.
+	(Routines should check for zeros first.)
+
+   FLT_UWORD_MIN
+	The bitmask of the smallest float above +0.  Call this number
+	REAL_FLT_MIN...
+
+   FLT_UWORD_EXP_MIN
+	The bitmask of the float representation of REAL_FLT_MIN's exponent.
+
+   FLT_UWORD_LOG_MIN
+	The bitmask of |log(REAL_FLT_MIN)|, rounding down.
+
+   FLT_SMALLEST_EXP
+	REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
+	-22 if they are).
+*/
+
+#ifdef _FLT_NO_DENORMALS
+#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
+#define FLT_UWORD_IS_SUBNORMAL(x) 0
+#define FLT_UWORD_MIN 0x00800000
+#define FLT_UWORD_EXP_MIN 0x42fc0000
+#define FLT_UWORD_LOG_MIN 0x42aeac50
+#define FLT_SMALLEST_EXP 1
+#else
+#define FLT_UWORD_IS_ZERO(x) ((x)==0)
+#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
+#define FLT_UWORD_MIN 0x00000001
+#define FLT_UWORD_EXP_MIN 0x43160000
+#define FLT_UWORD_LOG_MIN 0x42cff1b5
+#define FLT_SMALLEST_EXP -22
+#endif
+
 #ifdef __STDC__
 #define	__P(p)	p
 #else
 #define	__P(p)	()
 #endif
 
-#define	HUGE	((float)3.40282346638528860e+38)
-
 /* 
  * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
  * (one may replace the following line by "#include <values.h>")
diff --git a/newlib/libm/common/sf_cbrt.c b/newlib/libm/common/sf_cbrt.c
index c053d9548..fe632f0a8 100644
--- a/newlib/libm/common/sf_cbrt.c
+++ b/newlib/libm/common/sf_cbrt.c
@@ -53,13 +53,14 @@ G =  3.5714286566e-01; /* 5/14      = 0x3eb6db6e */
 	GET_FLOAT_WORD(hx,x);
 	sign=hx&0x80000000; 		/* sign= sign(x) */
 	hx  ^=sign;
-	if(hx>=0x7f800000) return(x+x); /* cbrt(NaN,INF) is itself */
-	if(hx==0) 
-	    return(x);		/* cbrt(0) is itself */
+	if(!FLT_UWORD_IS_FINITE(hx))
+	    return(x+x);		/* cbrt(NaN,INF) is itself */
+	if(FLT_UWORD_IS_ZERO(hx))
+	    return(x);			/* cbrt(0) is itself */
 
 	SET_FLOAT_WORD(x,hx);	/* x <- |x| */
     /* rough cbrt to 5 bits */
-	if(hx<0x00800000) 		/* subnormal number */
+	if(FLT_UWORD_IS_SUBNORMAL(hx)) 		/* subnormal number */
 	  {SET_FLOAT_WORD(t,0x4b800000); /* set t= 2**24 */
 	   t*=x; GET_FLOAT_WORD(high,t); SET_FLOAT_WORD(t,high/3+B2);
 	  }
diff --git a/newlib/libm/common/sf_expm1.c b/newlib/libm/common/sf_expm1.c
index e9108b7b8..4ba3b815a 100644
--- a/newlib/libm/common/sf_expm1.c
+++ b/newlib/libm/common/sf_expm1.c
@@ -27,7 +27,6 @@ static float
 one		= 1.0,
 huge		= 1.0e+30,
 tiny		= 1.0e-30,
-o_threshold	= 8.8721679688e+01,/* 0x42b17180 */
 ln2_hi		= 6.9313812256e-01,/* 0x3f317180 */
 ln2_lo		= 9.0580006145e-06,/* 0x3717f7d1 */
 invln2		= 1.4426950216e+00,/* 0x3fb8aa3b */
@@ -56,13 +55,12 @@ Q5  =  -2.0109921195e-07; /* 0xb457edbb */
 
     /* filter out huge and non-finite argument */
 	if(hx >= 0x4195b844) {			/* if |x|>=27*ln2 */
-	    if(hx >= 0x42b17218) {		/* if |x|>=88.721... */
-                if(hx>0x7f800000)
-		    return x+x; 	 /* NaN */
-		if(hx==0x7f800000)
-		    return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
-	        if(x > o_threshold) return huge*huge; /* overflow */
-	    }
+	    if(FLT_UWORD_IS_NAN(hx))
+	        return x+x;
+	    if(FLT_UWORD_IS_INFINITE(hx))
+		return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
+	    if(xsb == 0 && hx > FLT_UWORD_LOG_MAX) /* if x>=o_threshold */
+		return huge*huge; /* overflow */
 	    if(xsb!=0) { /* x < -27*ln2, return -1.0 with inexact */
 		if(x+tiny<(float)0.0)	/* raise inexact */
 		return tiny-one;	/* return -1 */
diff --git a/newlib/libm/common/sf_finite.c b/newlib/libm/common/sf_finite.c
index 4c48f400f..cb7026fbd 100644
--- a/newlib/libm/common/sf_finite.c
+++ b/newlib/libm/common/sf_finite.c
@@ -29,7 +29,8 @@
 {
 	__int32_t ix;
 	GET_FLOAT_WORD(ix,x);
-	return  (int)((__uint32_t)((ix&0x7fffffff)-0x7f800000)>>31);
+	ix &= 0x7fffffff;
+	return (FLT_UWORD_IS_FINITE(ix));
 }
 
 #ifdef _DOUBLE_IS_32BITS
diff --git a/newlib/libm/common/sf_ilogb.c b/newlib/libm/common/sf_ilogb.c
index ee65594b1..ec59406e1 100644
--- a/newlib/libm/common/sf_ilogb.c
+++ b/newlib/libm/common/sf_ilogb.c
@@ -27,15 +27,14 @@
 
 	GET_FLOAT_WORD(hx,x);
 	hx &= 0x7fffffff;
-	if(hx<0x00800000) {
-	    if(hx==0) 
-		return - INT_MAX;	/* ilogb(0) = 0x80000001 */
-	    else			/* subnormal x */
-	        for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
+	if(FLT_UWORD_IS_ZERO(hx))
+	    return - INT_MAX;	/* ilogb(0) = 0x80000001 */
+	if(FLT_UWORD_IS_SUBNORMAL(hx)) {
+	    for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
 	    return ix;
 	}
-	else if (hx<0x7f800000) return (hx>>23)-127;
-	else return INT_MAX;
+	else if (!FLT_UWORD_IS_FINITE(hx)) return INT_MAX;
+	else return (hx>>23)-127;
 }
 
 #ifdef _DOUBLE_IS_32BITS
diff --git a/newlib/libm/common/sf_log1p.c b/newlib/libm/common/sf_log1p.c
index 5ae7fb936..e09170f3e 100644
--- a/newlib/libm/common/sf_log1p.c
+++ b/newlib/libm/common/sf_log1p.c
@@ -51,6 +51,7 @@ static float zero = 0.0;
 	ax = hx&0x7fffffff;
 
 	k = 1;
+	if (!FLT_UWORD_IS_FINITE(hx)) return x+x;
 	if (hx < 0x3ed413d7) {			/* x < 0.41422  */
 	    if(ax>=0x3f800000) {		/* x <= -1.0 */
 		if(x==(float)-1.0) return -two25/zero; /* log1p(-1)=+inf */
@@ -65,8 +66,7 @@ static float zero = 0.0;
 	    }
 	    if(hx>0||hx<=((__int32_t)0xbe95f61f)) {
 		k=0;f=x;hu=1;}	/* -0.2929<x<0.41422 */
-	} 
-	if (hx >= 0x7f800000) return x+x;
+	}
 	if(k!=0) {
 	    if(hx<0x5a000000) {
 		u  = (float)1.0+x; 
diff --git a/newlib/libm/common/sf_logb.c b/newlib/libm/common/sf_logb.c
index 6e67637fe..f193f91f6 100644
--- a/newlib/libm/common/sf_logb.c
+++ b/newlib/libm/common/sf_logb.c
@@ -25,8 +25,8 @@
 	__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(FLT_UWORD_IS_ZERO(ix)) return (float)-1.0/fabsf(x);
+	if(!FLT_UWORD_IS_FINITE(ix)) return x*x;
 	if((ix>>=23)==0) 			/* IEEE 754 logb */
 		return -126.0; 
 	else
diff --git a/newlib/libm/common/sf_nan.c b/newlib/libm/common/sf_nan.c
index cb3e1cd03..c8d7027f8 100644
--- a/newlib/libm/common/sf_nan.c
+++ b/newlib/libm/common/sf_nan.c
@@ -21,3 +21,4 @@
 }
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
+
diff --git a/newlib/libm/common/sf_nextafter.c b/newlib/libm/common/sf_nextafter.c
index cd938d339..cea4da58d 100644
--- a/newlib/libm/common/sf_nextafter.c
+++ b/newlib/libm/common/sf_nextafter.c
@@ -29,15 +29,15 @@
 	ix = hx&0x7fffffff;		/* |x| */
 	iy = hy&0x7fffffff;		/* |y| */
 
-	if((ix>0x7f800000) ||   /* x is nan */ 
-	   (iy>0x7f800000))     /* y is nan */ 
-	   return x+y;				
+	if(FLT_UWORD_IS_NAN(ix) ||
+	   FLT_UWORD_IS_NAN(iy))
+	   return x+y;
 	if(x==y) return x;		/* x=y, return x */
-	if(ix==0) {				/* x == 0 */
-	    SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
+	if(FLT_UWORD_IS_ZERO(ix)) {		/* x == 0 */
+	    SET_FLOAT_WORD(x,(hy&0x80000000)|FLT_UWORD_MIN);
 	    y = x*x;
 	    if(y==x) return y; else return x;	/* raise underflow flag */
-	} 
+	}
 	if(hx>=0) {				/* x > 0 */
 	    if(hx>hy) {				/* x > y, x -= ulp */
 		hx -= 1;
@@ -52,7 +52,7 @@
 	    }
 	}
 	hy = hx&0x7f800000;
-	if(hy>=0x7f800000) return x+x;	/* overflow  */
+	if(hy>FLT_UWORD_MAX) return x+x;	/* overflow  */
 	if(hy<0x00800000) {		/* underflow */
 	    y = x*x;
 	    if(y!=x) {		/* raise underflow flag */
diff --git a/newlib/libm/common/sf_rint.c b/newlib/libm/common/sf_rint.c
index d38080a5d..6459b7a4c 100644
--- a/newlib/libm/common/sf_rint.c
+++ b/newlib/libm/common/sf_rint.c
@@ -33,15 +33,17 @@ TWO23[2]={
 #endif
 {
 	__int32_t i0,j0,sx;
-	__uint32_t i,i1;
+	__uint32_t i,i1,ix;
 	float t;
 	volatile float w;
 	GET_FLOAT_WORD(i0,x);
 	sx = (i0>>31)&1;
-	j0 = ((i0>>23)&0xff)-0x7f;
+	ix = (i0&0x7fffffff);
+	j0 = (ix>>23)-0x7f;
 	if(j0<23) {
-	    if(j0<0) { 	
-		if((i0&0x7fffffff)==0) return x;
+	    if(FLT_UWORD_IS_ZERO(ix))
+	        return x;
+	    if(j0<0) {
 		i1 = (i0&0x07fffff);
 		i0 &= 0xfff00000;
 		i0 |= ((i1|-i1)>>9)&0x400000;
@@ -58,8 +60,9 @@ TWO23[2]={
 		if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
 	    }
 	} else {
-	    if(j0==0x80) return x+x;	/* inf or NaN */
-	    else return x;		/* x is integral */
+	    if(!FLT_UWORD_IS_FINITE(ix)) return x+x; /* inf or NaN */
+	    else
+	      return x;		/* x is integral */
 	}
 	SET_FLOAT_WORD(x,i0);
 	w = TWO23[sx]+x;
diff --git a/newlib/libm/common/sf_scalbn.c b/newlib/libm/common/sf_scalbn.c
index fb67c7816..700060010 100644
--- a/newlib/libm/common/sf_scalbn.c
+++ b/newlib/libm/common/sf_scalbn.c
@@ -15,6 +15,7 @@
 
 #include "fdlibm.h"
 #include <limits.h>
+#include <float.h>
 
 #if INT_MAX > 50000
 #define OVERFLOW_INT 50000
@@ -40,25 +41,30 @@ tiny   = 1.0e-30;
 #endif
 {
 	__int32_t  k,ix;
+	__uint32_t hx;
+
 	GET_FLOAT_WORD(ix,x);
-        k = (ix&0x7f800000)>>23;		/* extract exponent */
-        if (k==0) {				/* 0 or subnormal x */
-            if ((ix&0x7fffffff)==0) return x; /* +-0 */
+	hx = ix&0x7fffffff;
+        k = hx>>23;		/* extract exponent */
+	if (FLT_UWORD_IS_ZERO(hx))
+	    return x;
+        if (!FLT_UWORD_IS_FINITE(hx))
+	    return x+x;		/* NaN or Inf */
+        if (FLT_UWORD_IS_SUBNORMAL(hx)) {
 	    x *= two25;
 	    GET_FLOAT_WORD(ix,x);
 	    k = ((ix&0x7f800000)>>23) - 25; 
             if (n< -50000) return tiny*x; 	/*underflow*/
-	    }
-        if (k==0xff) return x+x;		/* NaN or Inf */
+        }
         k = k+n; 
-        if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
+        if (k > FLT_LARGEST_EXP) return huge*copysignf(huge,x); /* overflow  */
         if (k > 0) 				/* normal result */
 	    {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
-        if (k <= -25) {
+        if (k < FLT_SMALLEST_EXP) {
             if (n > OVERFLOW_INT) 	/* in case integer overflow in n+k */
 		return huge*copysignf(huge,x);	/*overflow*/
 	    else return tiny*copysignf(tiny,x);	/*underflow*/
-       }
+        }
         k += 25;				/* subnormal result */
 	SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
         return x*twom25;