[PATCH] math: Optimize frexp with fast path for normal numbers
Adhemerval Zanella Netto
adhemerval.zanella@linaro.org
Tue Oct 21 17:27:29 GMT 2025
On 21/10/25 11:58, Osama Abdelkader wrote:
> 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>
Besides Paul's remarks, please split to one patch for each type (binary32/flt-32,
binary64/dbl-64, binary128/ldbl-128, and intel96/ldbl-96).
> ---
> 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);
Could you also use the macros from math_config.h to obtain the bit representation
from the floating-point number? It can be written as:
int64_t ix = asuint64 (x);
> - int32_t ex = 0x7ff & (ix >> 52);
> - int e = 0;
> + uint32_t ex = 0x7ff & (ix >> 52);
It can be either:
uint32_t ex = get_exponent (ix) + EXPONENT_BIAS;
Or maybe add a new function get_exponent_unbiased.
>
> - 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));
I think this like is too long, maybe:
return asdouble ((ix & (MANTISSA_MASK | SIGN_MASK)
| INT64_C (0x3fe0000000000000));
And maybe add INT64_C (0x3fe0000000000000) as a macro.
> + }
>
> + /* 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);
ix = asuint64 (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;
More information about the Libc-alpha
mailing list