[PATCH v2 1/4] math: Remove UB and optimize double ilogb
Adhemerval Zanella
adhemerval.zanella@linaro.org
Mon Apr 28 17:03:41 GMT 2025
The subnormal exponent calculation invokes UB by left shifting the
signed expoenent to find the first leading bit. The implementation
also used 32 bits operations, which is generates suboptimal code
in 64 bits architectures.
The patch reimplements ilogb using the math_config.h macros and
uses the new stdbit function to simplify the subnormal handling.
On aarch64 it generates better code:
* master:
0000000000000000 <__ieee754_ilogb>:
0: 9e660000 fmov x0, d0
4: d360fc02 lsr x2, x0, #32
8: d360f801 ubfx x1, x0, #32, #31
c: f26c285f tst x2, #0x7ff00000
10: 540001a1 b.ne 44 <__ieee754_ilogb+0x44> // b.any
14: 2a000022 orr w2, w1, w0
18: 34000322 cbz w2, 7c <__ieee754_ilogb+0x7c>
1c: 35000221 cbnz w1, 60 <__ieee754_ilogb+0x60>
20: 2a0003e1 mov w1, w0
24: 7100001f cmp w0, #0x0
28: 12808240 mov w0, #0xfffffbed // #-1043
2c: 540000ad b.le 40 <__ieee754_ilogb+0x40>
30: 531f7821 lsl w1, w1, #1
34: 51000400 sub w0, w0, #0x1
38: 7100003f cmp w1, #0x0
3c: 54ffffac b.gt 30 <__ieee754_ilogb+0x30>
40: d65f03c0 ret
44: 13147c20 asr w0, w1, #20
48: 12b00202 mov w2, #0x7fefffff // #2146435071
4c: 510ffc00 sub w0, w0, #0x3ff
50: 6b02003f cmp w1, w2
54: 12b00001 mov w1, #0x7fffffff // #2147483647
58: 1a819000 csel w0, w0, w1, ls // ls = plast
5c: d65f03c0 ret
60: 53155021 lsl w1, w1, #11
64: 12807fa0 mov w0, #0xfffffc02 // #-1022
68: 531f7821 lsl w1, w1, #1
6c: 51000400 sub w0, w0, #0x1
70: 7100003f cmp w1, #0x0
74: 54ffffac b.gt 68 <__ieee754_ilogb+0x68>
78: d65f03c0 ret
7c: 320107e0 mov w0, #0x80000001 // #-2147483647
80: d65f03c0 ret
* patch:
0000000000000000 <__ieee754_ilogb>:
0: 9e660001 fmov x1, d0
4: d374f820 ubfx x0, x1, #52, #11
8: 350000e0 cbnz w0, 24 <__ieee754_ilogb+0x24>
c: d374cc21 lsl x1, x1, #12
10: b4000141 cbz x1, 38 <__ieee754_ilogb+0x38>
14: dac01021 clz x1, x1
18: 12807fc0 mov w0, #0xfffffc01 // #-1023
1c: 4b010000 sub w0, w0, w1
20: d65f03c0 ret
24: 711ffc1f cmp w0, #0x7ff
28: 510ffc00 sub w0, w0, #0x3ff
2c: 12b00001 mov w1, #0x7fffffff // #2147483647
30: 1a811000 csel w0, w0, w1, ne // ne = any
34: d65f03c0 ret
38: 320107e0 mov w0, #0x80000001 // #-2147483647
3c: d65f03c0 ret
Other architecture with support for stdc_leading_zeros and/or
__builtin_clzll should have similar improvements.
Checked on aarch64-linux-gnu and x86_64-linux-gnu.
---
sysdeps/ieee754/dbl-64/e_ilogb.c | 80 ++++++++++++--------------------
1 file changed, 29 insertions(+), 51 deletions(-)
diff --git a/sysdeps/ieee754/dbl-64/e_ilogb.c b/sysdeps/ieee754/dbl-64/e_ilogb.c
index 1e338a59c1..89e7498266 100644
--- a/sysdeps/ieee754/dbl-64/e_ilogb.c
+++ b/sysdeps/ieee754/dbl-64/e_ilogb.c
@@ -1,63 +1,41 @@
-/* @(#)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.
- * ====================================================
- */
+/* Get integer exponent of a floating-point value.
+ Copyright (C) 1999-2025 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
-#endif
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
-/* ilogb(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)
- */
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <https://www.gnu.org/licenses/>. */
#include <limits.h>
#include <math.h>
-#include <math_private.h>
+#include <stdbit.h>
+#include "math_config.h"
int
__ieee754_ilogb (double x)
{
- int32_t hx, lx, ix;
-
- GET_HIGH_WORD (hx, x);
- hx &= 0x7fffffff;
- if (hx < 0x00100000)
+ uint64_t ux = asuint64 (x);
+ int ex = (ux & ~SIGN_MASK) >> MANTISSA_WIDTH;
+ if (ex == 0) /* zero or subnormal */
{
- GET_LOW_WORD (lx, x);
- if ((hx | lx) == 0)
- return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
- else /* subnormal x */
- if (hx == 0)
- {
- for (ix = -1043; lx > 0; lx <<= 1)
- ix -= 1;
- }
- else
- {
- for (ix = -1022, hx <<= 11; hx > 0; hx <<= 1)
- ix -= 1;
- }
- return ix;
+ /* Clear sign and exponent */
+ ux <<= 12;
+ if (ux == 0)
+ return FP_ILOGB0;
+ /* subnormal */
+ return -1023 - stdc_leading_zeros (ux);
}
- else if (hx < 0x7ff00000)
- return (hx >> 20) - 1023;
- else if (FP_ILOGBNAN != INT_MAX)
- {
- /* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
- GET_LOW_WORD (lx, x);
- if (((hx ^ 0x7ff00000) | lx) == 0)
- return INT_MAX;
- }
- return FP_ILOGBNAN;
+ if (ex == EXPONENT_MASK >> MANTISSA_WIDTH) /* NaN or Inf */
+ return ux << 12 ? FP_ILOGBNAN : INT_MAX;
+ return ex - 1023;
}
--
2.43.0
More information about the Libc-alpha
mailing list