[PATCH] math: Optimize frexp with fast path for normal numbers
Osama Abdelkader
osama.abdelkader@gmail.com
Tue Oct 21 14:58:17 GMT 2025
Add fast path optimization to frexp family functions for normal floating-point
numbers, which represent most of real-world usage. The optimization uses a
single unsigned comparison to identify normal numbers and returns immediately
via bit manipulation, avoiding branches and floating-point operations.
Performance improvements:
- Fewer branches
- Better branch prediction: __glibc_likely hints for common case
The key insight is the check (ex - 1) < MAX_NORMAL_EXPONENT with unsigned
arithmetic, which elegantly identifies normal numbers while excluding zero,
subnormal, infinity, and NaN values.
Files modified:
- sysdeps/ieee754/dbl-64/s_frexp.c
- sysdeps/ieee754/flt-32/s_frexpf.c
- sysdeps/ieee754/ldbl-96/s_frexpl.c
- sysdeps/ieee754/ldbl-128/s_frexpl.c
All existing tests pass, maintaining exact correctness for all edge cases.
Suggested-by: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>
---
sysdeps/ieee754/dbl-64/s_frexp.c | 29 ++++++++++++++------------
sysdeps/ieee754/flt-32/s_frexpf.c | 31 ++++++++++++++++++----------
sysdeps/ieee754/ldbl-128/s_frexpl.c | 32 +++++++++++++++++++----------
sysdeps/ieee754/ldbl-96/s_frexpl.c | 29 +++++++++++++++++---------
4 files changed, 76 insertions(+), 45 deletions(-)
diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c
index 9c45819d4d..d873790f32 100644
--- a/sysdeps/ieee754/dbl-64/s_frexp.c
+++ b/sysdeps/ieee754/dbl-64/s_frexp.c
@@ -18,6 +18,7 @@
#include <inttypes.h>
#include <math.h>
#include <math_private.h>
+#include "math_config.h"
#include <libm-alias-double.h>
/*
@@ -36,22 +37,24 @@ __frexp (double x, int *eptr)
{
int64_t ix;
EXTRACT_WORDS64 (ix, x);
- int32_t ex = 0x7ff & (ix >> 52);
- int e = 0;
+ uint32_t ex = 0x7ff & (ix >> 52);
- if (__glibc_likely (ex != 0x7ff && x != 0.0))
+ /* Fast path for normal numbers. */
+ if (__glibc_likely ((ex - 1) < 0x7fe))
{
- /* Not zero and finite. */
- e = ex - 1022;
- if (__glibc_unlikely (ex == 0))
- {
- /* Subnormal. */
- x *= 0x1p54;
- EXTRACT_WORDS64 (ix, x);
- ex = 0x7ff & (ix >> 52);
- e = ex - 1022 - 54;
- }
+ *eptr = ex - 1022;
+ return asdouble ((ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000));
+ }
+ /* Handle special cases: zero, subnormal, infinity, NaN. */
+ int e = 0;
+ if (__glibc_likely (ex != 0x7ff && x != 0.0))
+ {
+ /* Subnormal. */
+ x *= 0x1p54;
+ EXTRACT_WORDS64 (ix, x);
+ ex = 0x7ff & (ix >> 52);
+ e = ex - 1022 - 54;
ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
INSERT_WORDS64 (x, ix);
}
diff --git a/sysdeps/ieee754/flt-32/s_frexpf.c b/sysdeps/ieee754/flt-32/s_frexpf.c
index 59fef66a6b..46f4cb2497 100644
--- a/sysdeps/ieee754/flt-32/s_frexpf.c
+++ b/sysdeps/ieee754/flt-32/s_frexpf.c
@@ -18,6 +18,7 @@ static char rcsid[] = "$NetBSD: s_frexpf.c,v 1.5 1995/05/10 20:47:26 jtc Exp $";
#include <math.h>
#include <math_private.h>
+#include "math_config.h"
#include <libm-alias-float.h>
static const float
@@ -25,19 +26,27 @@ two25 = 3.3554432000e+07; /* 0x4c000000 */
float __frexpf(float x, int *eptr)
{
- int32_t hx,ix;
+ int32_t hx;
GET_FLOAT_WORD(hx,x);
- ix = 0x7fffffff&hx;
+ uint32_t ex = 0xff & (hx >> 23);
+
+ /* Fast path for normal numbers. */
+ if (__glibc_likely ((ex - 1) < 0xfe))
+ {
+ *eptr = ex - 126;
+ return asfloat ((hx & 0x807fffff) | 0x3f000000);
+ }
+
+ /* Handle special cases: zero, subnormal, infinity, NaN. */
+ int32_t ix = 0x7fffffff & hx;
*eptr = 0;
- if(ix>=0x7f800000||(ix==0)) return x + x; /* 0,inf,nan */
- if (ix<0x00800000) { /* subnormal */
- x *= two25;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff;
- *eptr = -25;
- }
- *eptr += (ix>>23)-126;
- hx = (hx&0x807fffff)|0x3f000000;
+ if (ix >= 0x7f800000 || (ix == 0)) return x + x; /* 0,inf,nan */
+ /* Subnormal */
+ x *= two25;
+ GET_FLOAT_WORD(hx,x);
+ ix = hx & 0x7fffffff;
+ *eptr = (ix >> 23) - 126 - 25;
+ hx = (hx & 0x807fffff) | 0x3f000000;
SET_FLOAT_WORD(x,hx);
return x;
}
diff --git a/sysdeps/ieee754/ldbl-128/s_frexpl.c b/sysdeps/ieee754/ldbl-128/s_frexpl.c
index e4db093a05..1649a5a20a 100644
--- a/sysdeps/ieee754/ldbl-128/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-128/s_frexpl.c
@@ -35,19 +35,29 @@ two114 = L(2.0769187434139310514121985316880384E+34); /* 0x4071000000000000, 0 *
_Float128 __frexpl(_Float128 x, int *eptr)
{
- uint64_t hx, lx, ix;
+ uint64_t hx, lx;
GET_LDOUBLE_WORDS64(hx,lx,x);
- ix = 0x7fffffffffffffffULL&hx;
+ uint64_t ex = 0x7fff & (hx >> 48);
+
+ /* Fast path for normal numbers. */
+ if (__glibc_likely ((ex - 1U) < 0x7ffe))
+ {
+ *eptr = ex - 16382;
+ hx = (hx & 0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
+ SET_LDOUBLE_MSW64(x,hx);
+ return x;
+ }
+
+ /* Handle special cases: zero, subnormal, infinity, NaN. */
+ uint64_t ix = 0x7fffffffffffffffULL & hx;
*eptr = 0;
- if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x + x;/* 0,inf,nan */
- if (ix<0x0001000000000000ULL) { /* subnormal */
- x *= two114;
- GET_LDOUBLE_MSW64(hx,x);
- ix = hx&0x7fffffffffffffffULL;
- *eptr = -114;
- }
- *eptr += (ix>>48)-16382;
- hx = (hx&0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
+ if (ix >= 0x7fff000000000000ULL || ((ix | lx) == 0)) return x + x;/* 0,inf,nan */
+ /* Subnormal */
+ x *= two114;
+ GET_LDOUBLE_MSW64(hx,x);
+ ix = hx & 0x7fffffffffffffffULL;
+ *eptr = (ix >> 48) - 16382 - 114;
+ hx = (hx & 0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
SET_LDOUBLE_MSW64(x,hx);
return x;
}
diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c
index c610704cca..1e99c6ed9c 100644
--- a/sysdeps/ieee754/ldbl-96/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c
@@ -41,18 +41,27 @@ two65 = 3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */
long double __frexpl(long double x, int *eptr)
{
- uint32_t se, hx, ix, lx;
+ uint32_t se, hx, lx;
GET_LDOUBLE_WORDS(se,hx,lx,x);
- ix = 0x7fff&se;
+ uint32_t ex = 0x7fff & se;
+
+ /* Fast path for normal numbers. */
+ if (__glibc_likely ((ex - 1U) < 0x7ffe))
+ {
+ *eptr = ex - 16382;
+ se = (se & 0x8000) | 0x3ffe;
+ SET_LDOUBLE_EXP(x,se);
+ return x;
+ }
+
+ /* Handle special cases: zero, subnormal, infinity, NaN. */
*eptr = 0;
- if(ix==0x7fff||((ix|hx|lx)==0)) return x + x; /* 0,inf,nan */
- if (ix==0x0000) { /* subnormal */
- x *= two65;
- GET_LDOUBLE_EXP(se,x);
- ix = se&0x7fff;
- *eptr = -65;
- }
- *eptr += ix-16382;
+ if (ex == 0x7fff || ((ex | hx | lx) == 0)) return x + x; /* 0,inf,nan */
+ /* Subnormal */
+ x *= two65;
+ GET_LDOUBLE_EXP(se,x);
+ ex = se & 0x7fff;
+ *eptr = ex - 16382 - 65;
se = (se & 0x8000) | 0x3ffe;
SET_LDOUBLE_EXP(x,se);
return x;
--
2.43.0
More information about the Libc-alpha
mailing list