Fix ilogb(3) implementation. The results have to match FP_ILOGB0 and

FP_ILOGBNAN which isn't the case for the amd64 and i386 assembly versions.
Drop these in favour of C implementations.  Als reimplement ilogbl(3)
by providing separate ld80 and ld128 implementations that replace the
existing implementation which may hit an infinite loop when built for
quad-precision long double.

ok patrick@, gkoehler@

5 files changed, 110 insertions, 85 deletions

diff --git a/lib/libm/src/ld128/s_ilogbl.c b/lib/libm/src/ld128/s_ilogbl.c
new file mode 100644
index 00000000000..879778e2f7d
--- /dev/null
+++ b/lib/libm/src/ld128/s_ilogbl.c
@@ -0,0 +1,51 @@
+/* s_ilogbl.c -- long double version of s_ilogb.c.
+ * Conversion to 128-bit long double by Mark Kettenis, kettenis@openbsd.org.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* ilogbl(long double x)
+ * return the binary exponent of non-zero x
+ * ilogb(0) = FP_ILOGB0
+ * ilogb(NaN) = FP_ILOGBNAN (no signal is raised)
+ * ilogb(inf) = INT_MAX (no signal is raised)
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+int
+ilogbl(long double x)
+{
+	int64_t hx,lx;
+	int32_t ix;
+
+	GET_LDOUBLE_WORDS64(hx,lx,x);
+	hx &= 0x7fffffffffffffffLL;
+	if(hx<0x0001000000000000LL) {
+	    if((hx|lx)==0) 
+		return FP_ILOGB0;	/* ilogb(0) = FP_ILOGB0 */
+	    else			/* subnormal x */
+		if(hx==0) {
+		    for (ix = -16431; lx>0; lx<<=1) ix -=1;
+		} else {
+		    for (ix = -16382, hx<<=15; hx>0; hx<<=1) ix -=1;
+		}
+	    return ix;
+	}
+	else if (hx<0x7fff000000000000LL) return (hx>>48)-16383;
+	else if (hx>0x7fff000000000000LL || lx!=0) return FP_ILOGBNAN;
+	else return INT_MAX;
+}
+DEF_STD(ilogbl);
diff --git a/lib/libm/src/ld80/s_ilogbl.c b/lib/libm/src/ld80/s_ilogbl.c
new file mode 100644
index 00000000000..dba17022009
--- /dev/null
+++ b/lib/libm/src/ld80/s_ilogbl.c
@@ -0,0 +1,50 @@
+/* s_ilogbl.c -- long double version of s_ilogb.c.
+ * Conversion to 80-bit long double by Mark Kettenis, kettenis@openbsd.org.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* ilogbl(long double x)
+ * return the binary exponent of non-zero x
+ * ilogb(0) = FP_ILOGB0
+ * ilogb(NaN) = FP_ILOGBNAN (no signal is raised)
+ * ilogb(inf) = INT_MAX (no signal is raised)
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+int
+ilogbl(long double x)
+{
+	int32_t esx,hx,lx,ix;
+
+	GET_LDOUBLE_WORDS(esx,hx,lx,x);
+	esx &= 0x7fff;
+	if(esx==0) {
+	    if((hx|lx)==0) 
+		return FP_ILOGB0;	/* ilogb(0) = FP_ILOGB0 */
+	    else			/* subnormal x */
+		if(hx==0) {
+		    for (ix = -16414; lx>0; lx<<=1) ix -=1;
+		} else {
+		    for (ix = -16382; hx>0; hx<<=1) ix -=1;
+		}
+	    return ix;
+	}
+	else if (esx<0x7fff) return (esx)-16383;
+	else if ((hx&0x7fffffff|lx)!=0) return FP_ILOGBNAN;
+	else return INT_MAX;
+}
+DEF_STD(ilogbl);
diff --git a/lib/libm/src/s_ilogb.c b/lib/libm/src/s_ilogb.c
index 14dec8dacb5..bfa627eb7e5 100644
--- a/lib/libm/src/s_ilogb.c
+++ b/lib/libm/src/s_ilogb.c
@@ -12,8 +12,9 @@
 
 /* ilogb(double x)
  * return the binary exponent of non-zero x
- * ilogb(0) = 0x80000001
- * ilogb(inf/NaN) = 0x7fffffff (no signal is raised)
+ * ilogb(0) = FP_ILOGB0
+ * ilogb(NaN) = FP_ILOGBNAN (no signal is raised)
+ * ilogb(inf) = INT_MAX (no signal is raised)
  */
 
 #include <float.h>
@@ -31,7 +32,7 @@ ilogb(double x)
 	if(hx<0x00100000) {
 	    GET_LOW_WORD(lx,x);
 	    if((hx|lx)==0) 
-		return 0x80000001;	/* ilogb(0) = 0x80000001 */
+		return FP_ILOGB0;	/* ilogb(0) = FP_ILOGB0 */
 	    else			/* subnormal x */
 		if(hx==0) {
 		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
@@ -41,7 +42,8 @@ ilogb(double x)
 	    return ix;
 	}
 	else if (hx<0x7ff00000) return (hx>>20)-1023;
-	else return 0x7fffffff;
+	else if (hx>0x7ff00000 || lx!=0) return FP_ILOGBNAN;
+	else return INT_MAX;
 }
 DEF_STD(ilogb);
 LDBL_MAYBE_CLONE(ilogb);
diff --git a/lib/libm/src/s_ilogbf.c b/lib/libm/src/s_ilogbf.c
index 04c3e316d54..5f2816382d5 100644
--- a/lib/libm/src/s_ilogbf.c
+++ b/lib/libm/src/s_ilogbf.c
@@ -25,12 +25,13 @@ ilogbf(float x)
 	hx &= 0x7fffffff;
 	if(hx<0x00800000) {
 	    if(hx==0) 
-		return 0x80000001;	/* ilogb(0) = 0x80000001 */
+		return FP_ILOGB0;	/* ilogb(0) = FP_ILOGB */
 	    else			/* subnormal x */
 	        for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
 	    return ix;
 	}
 	else if (hx<0x7f800000) return (hx>>23)-127;
-	else return 0x7fffffff;
+	else if (hx>0x7f800000) return FP_ILOGBNAN;
+	else return INT_MAX;
 }
 DEF_STD(ilogbf);
diff --git a/lib/libm/src/s_ilogbl.c b/lib/libm/src/s_ilogbl.c
deleted file mode 100644
index 9c1f83fff4f..00000000000
--- a/lib/libm/src/s_ilogbl.c
+++ /dev/null
@@ -1,79 +0,0 @@
-/*	$OpenBSD: s_ilogbl.c,v 1.2 2016/09/12 19:47:02 guenther Exp $	*/
-/*
- * From: @(#)s_ilogb.c 5.1 93/09/24
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#include <sys/types.h>
-#include <machine/ieee.h>
-#include <float.h>
-#include <limits.h>
-#include <math.h>
-
-int
-ilogbl(long double x)
-{
-	struct ieee_ext *p = (struct ieee_ext *)&x;
-	unsigned long m;
-	int b;
-
-	if (p->ext_exp == 0) {
-		if ((p->ext_fracl
-#ifdef EXT_FRACLMBITS
-			| p->ext_fraclm
-#endif /* EXT_FRACLMBITS */
-#ifdef EXT_FRACHMBITS
-			| p->ext_frachm
-#endif /* EXT_FRACHMBITS */
-			| p->ext_frach) == 0)
-			return (FP_ILOGB0);
-		/* denormalized */
-		if (p->ext_frach == 0
-#ifdef EXT_FRACHMBITS
-			&& p->ext_frachm == 0
-#endif
-			) {
-			m = 1lu << (EXT_FRACLBITS - 1);
-			for (b = EXT_FRACHBITS; !(p->ext_fracl & m); m >>= 1)
-				b++;
-#if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS)
-			m = 1lu << (EXT_FRACLMBITS - 1);
-			for (b += EXT_FRACHMBITS; !(p->ext_fraclm & m); m >>= 1)
-				b++;
-#endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */
-		} else {
-			m = 1lu << (EXT_FRACHBITS - 1);
-			for (b = 0; !(p->ext_frach & m); m >>= 1)
-				b++;
-#ifdef EXT_FRACHMBITS
-			m = 1lu << (EXT_FRACHMBITS - 1);
-			for (; !(p->ext_frachm & m); m >>= 1)
-				b++;
-#endif /* EXT_FRACHMBITS */
-		}
-#ifdef EXT_IMPLICIT_NBIT
-		b++;
-#endif
-		return (LDBL_MIN_EXP - b - 1);
-	} else if (p->ext_exp < (LDBL_MAX_EXP << 1) - 1)
-		return (p->ext_exp - LDBL_MAX_EXP + 1);
-	else if (p->ext_fracl != 0
-#ifdef EXT_FRACLMBITS
-		|| p->ext_fraclm != 0
-#endif /* EXT_FRACLMBITS */
-#ifdef EXT_FRACHMBITS
-		|| p->ext_frachm != 0
-#endif /* EXT_FRACHMBITS */
-		|| p->ext_frach != 0)
-		return (FP_ILOGBNAN);
-	else
-		return (INT_MAX);
-}
-DEF_STD(ilogbl);