As reported on bugs@ our ilogb(3) implementation is somewhat buggy.
There are several issues:
- The amd64 and i386 assembler versions not only return results that
don't match the FP_ILOGB0 and FP_ILOGBNAN defines in <math.h> but
ultimately return results that are incompatible with C99 and C11.
- The C versions don't use the FP_ILOGB0 and FP_ILOGBNAN defines.
- The ilogbl(3) implementation turns into an infinite loop with
182-bit long doubles in some cases.
Here is a diff that fixes those issues by:
- Dropping the amd64 and i386 versions. Fixing the corner cases in
assembly is hard, and the C implementation should be fast enough for
regular floating-point values.
- Replaces the broken C implementation of ilogbl(3) with separate
128-bit and 80-bit long double implementations.
- Add the regress test from FreeBSD.
Technically, this should be a libm major bump since the behaviour
changes on amd64 and i386. But given that the corner cases were so
blatantly wrong on those platforms, I don't think that's necessary.
ok?
Index: lib/libm/Makefile
===================================================================
RCS file: /cvs/src/lib/libm/Makefile,v
retrieving revision 1.120
diff -u -p -r1.120 Makefile
--- lib/libm/Makefile 28 Jun 2020 08:22:57 -0000 1.120
+++ lib/libm/Makefile 31 Oct 2020 14:52:14 -0000
@@ -23,7 +23,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S
invtrig.c \
s_atan.S s_atanf.S s_ceil.S s_ceilf.S s_copysign.S s_copysignf.S \
s_cos.S s_cosf.S s_floor.S s_floorf.S \
- s_ilogb.S s_ilogbf.S s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
+ s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S s_rint.S s_rintf.S\
s_scalbnf.S s_significand.S s_significandf.S \
s_sin.S s_sinf.S s_tan.S s_tanf.S
@@ -36,7 +36,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S
invtrig.c \
s_atan.S s_atanf.S s_ceil.S s_ceilf.S s_copysign.S s_copysignf.S \
s_cos.S s_cosf.S s_floor.S s_floorf.S \
- s_ilogb.S s_ilogbf.S s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
+ s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S \
s_rint.S s_rintf.S s_scalbnf.S s_significand.S \
s_significandf.S s_sin.S s_sinf.S s_tan.S s_tanf.S
Index: lib/libm/arch/amd64/s_ilogb.S
===================================================================
RCS file: lib/libm/arch/amd64/s_ilogb.S
diff -N lib/libm/arch/amd64/s_ilogb.S
--- lib/libm/arch/amd64/s_ilogb.S 3 Jul 2018 22:43:34 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,20 +0,0 @@
-/* $OpenBSD: s_ilogb.S,v 1.5 2018/07/03 22:43:34 mortimer Exp $ */
-/*
- * Written by J.T. Conklin <[email protected]>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include "abi.h"
-
-ENTRY(ilogb)
- RETGUARD_SETUP(ilogb, r11)
- movsd %xmm0,-8(%rsp)
- fldl -8(%rsp)
- fxtract
- fstp %st
- fistpl -8(%rsp)
- movl -8(%rsp),%eax
- RETGUARD_CHECK(ilogb, r11)
- ret
-END_STD(ilogb)
Index: lib/libm/arch/amd64/s_ilogbf.S
===================================================================
RCS file: lib/libm/arch/amd64/s_ilogbf.S
diff -N lib/libm/arch/amd64/s_ilogbf.S
--- lib/libm/arch/amd64/s_ilogbf.S 3 Jul 2018 22:43:34 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,20 +0,0 @@
-/* $OpenBSD: s_ilogbf.S,v 1.5 2018/07/03 22:43:34 mortimer Exp $ */
-/*
- * Written by J.T. Conklin <[email protected]>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include "abi.h"
-
-ENTRY(ilogbf)
- RETGUARD_SETUP(ilogbf, r11)
- movss %xmm0,-4(%rsp)
- flds -4(%rsp)
- fxtract
- fstp %st
- fistpl -4(%rsp)
- movl -4(%rsp),%eax
- RETGUARD_CHECK(ilogbf, r11)
- ret
-END_STD(ilogbf)
Index: lib/libm/arch/i387/s_ilogb.S
===================================================================
RCS file: lib/libm/arch/i387/s_ilogb.S
diff -N lib/libm/arch/i387/s_ilogb.S
--- lib/libm/arch/i387/s_ilogb.S 12 Sep 2016 19:47:02 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,23 +0,0 @@
-/* $OpenBSD: s_ilogb.S,v 1.6 2016/09/12 19:47:02 guenther Exp $ */
-/*
- * Written by J.T. Conklin <[email protected]>.
- * Public domain.
- */
-
-#include "DEFS.h"
-
-ENTRY(ilogb)
- pushl %ebp
- movl %esp,%ebp
- subl $4,%esp
-
- fldl 8(%ebp)
- fxtract
- fstp %st
-
- fistpl -4(%ebp)
- movl -4(%ebp),%eax
-
- leave
- ret
-END_STD(ilogb)
Index: lib/libm/arch/i387/s_ilogbf.S
===================================================================
RCS file: lib/libm/arch/i387/s_ilogbf.S
diff -N lib/libm/arch/i387/s_ilogbf.S
--- lib/libm/arch/i387/s_ilogbf.S 12 Sep 2016 19:47:02 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,23 +0,0 @@
-/* $OpenBSD: s_ilogbf.S,v 1.6 2016/09/12 19:47:02 guenther Exp $ */
-/*
- * Written by J.T. Conklin <[email protected]>.
- * Public domain.
- */
-
-#include "DEFS.h"
-
-ENTRY(ilogbf)
- pushl %ebp
- movl %esp,%ebp
- subl $4,%esp
-
- flds 8(%ebp)
- fxtract
- fstp %st
-
- fistpl -4(%ebp)
- movl -4(%ebp),%eax
-
- leave
- ret
-END_STD(ilogbf)
Index: lib/libm/src/s_ilogb.c
===================================================================
RCS file: /cvs/src/lib/libm/src/s_ilogb.c,v
retrieving revision 1.11
diff -u -p -r1.11 s_ilogb.c
--- lib/libm/src/s_ilogb.c 12 Sep 2016 19:47:02 -0000 1.11
+++ lib/libm/src/s_ilogb.c 31 Oct 2020 14:52:14 -0000
@@ -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);
Index: lib/libm/src/s_ilogbf.c
===================================================================
RCS file: /cvs/src/lib/libm/src/s_ilogbf.c,v
retrieving revision 1.4
diff -u -p -r1.4 s_ilogbf.c
--- lib/libm/src/s_ilogbf.c 12 Sep 2016 19:47:02 -0000 1.4
+++ lib/libm/src/s_ilogbf.c 31 Oct 2020 14:52:14 -0000
@@ -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);
Index: lib/libm/src/s_ilogbl.c
===================================================================
RCS file: lib/libm/src/s_ilogbl.c
diff -N lib/libm/src/s_ilogbl.c
--- lib/libm/src/s_ilogbl.c 12 Sep 2016 19:47:02 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -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);
Index: lib/libm/src/ld128/s_ilogbl.c
===================================================================
RCS file: lib/libm/src/ld128/s_ilogbl.c
diff -N lib/libm/src/ld128/s_ilogbl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/ld128/s_ilogbl.c 31 Oct 2020 14:52:14 -0000
@@ -0,0 +1,51 @@
+/* s_ilogbl.c -- long double version of s_ilogb.c.
+ * Conversion to 128-bit long double by Mark Kettenis, [email protected].
+ */
+
+/*
+ * ====================================================
+ * 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);
Index: lib/libm/src/ld80/s_ilogbl.c
===================================================================
RCS file: lib/libm/src/ld80/s_ilogbl.c
diff -N lib/libm/src/ld80/s_ilogbl.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ lib/libm/src/ld80/s_ilogbl.c 31 Oct 2020 14:52:14 -0000
@@ -0,0 +1,50 @@
+/* s_ilogbl.c -- long double version of s_ilogb.c.
+ * Conversion to 80-bit long double by Mark Kettenis, [email protected].
+ */
+
+/*
+ * ====================================================
+ * 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);
Index: regress/lib/libm/msun/Makefile
===================================================================
RCS file: /cvs/src/regress/lib/libm/msun/Makefile,v
retrieving revision 1.1.1.1
diff -u -p -r1.1.1.1 Makefile
--- regress/lib/libm/msun/Makefile 21 Feb 2019 16:14:03 -0000 1.1.1.1
+++ regress/lib/libm/msun/Makefile 31 Oct 2020 14:52:14 -0000
@@ -3,6 +3,7 @@
TESTS =
TESTS += conj_test
TESTS += fenv_test
+TESTS += ilogb_test
TESTS += lrint_test
PROGS= ${TESTS}
Index: regress/lib/libm/msun/ilogb_test.c
===================================================================
RCS file: regress/lib/libm/msun/ilogb_test.c
diff -N regress/lib/libm/msun/ilogb_test.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regress/lib/libm/msun/ilogb_test.c 31 Oct 2020 14:52:14 -0000
@@ -0,0 +1,83 @@
+/*-
+ * Copyright (c) 2004 Stefan Farfeleder
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: head/lib/msun/tests/ilogb_test.c 292328 2015-12-16 09:11:11Z ngie
$
+ */
+
+#include <assert.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+int
+main(void)
+{
+ char buf[128], *end;
+ double d;
+ float f;
+ long double ld;
+ int e, i;
+
+ printf("1..3\n");
+ assert(ilogb(0) == FP_ILOGB0);
+ assert(ilogb(NAN) == FP_ILOGBNAN);
+ assert(ilogb(INFINITY) == INT_MAX);
+ for (e = DBL_MIN_EXP - DBL_MANT_DIG; e < DBL_MAX_EXP; e++) {
+ snprintf(buf, sizeof(buf), "0x1.p%d", e);
+ d = strtod(buf, &end);
+ assert(*end == '\0');
+ i = ilogb(d);
+ assert(i == e);
+ }
+ printf("ok 1 - ilogb\n");
+
+ assert(ilogbf(0) == FP_ILOGB0);
+ assert(ilogbf(NAN) == FP_ILOGBNAN);
+ assert(ilogbf(INFINITY) == INT_MAX);
+ for (e = FLT_MIN_EXP - FLT_MANT_DIG; e < FLT_MAX_EXP; e++) {
+ snprintf(buf, sizeof(buf), "0x1.p%d", e);
+ f = strtof(buf, &end);
+ assert(*end == '\0');
+ i = ilogbf(f);
+ assert(i == e);
+ }
+ printf("ok 2 - ilogbf\n");
+
+ assert(ilogbl(0) == FP_ILOGB0);
+ assert(ilogbl(NAN) == FP_ILOGBNAN);
+ assert(ilogbl(INFINITY) == INT_MAX);
+ for (e = LDBL_MIN_EXP - LDBL_MANT_DIG; e < LDBL_MAX_EXP; e++) {
+ snprintf(buf, sizeof(buf), "0x1.p%d", e);
+ ld = strtold(buf, &end);
+ assert(*end == '\0');
+ i = ilogbl(ld);
+ assert(i == e);
+ }
+ printf("ok 3 - ilogbl\n");
+
+ return (0);
+}
