From e05476b5c8ef33cf4a10332663ad06a10507bc2b Mon Sep 17 00:00:00 2001 From: Osama Abdelkader Date: Thu, 23 Oct 2025 18:06:29 +0300 Subject: [PATCH] math: Optimize frexp (binary64) with fast path for normal numbers Add fast path optimization for frexp using a single unsigned comparison to identify normal floating-point numbers and return immediately via arithmetic on the bit representation. The implementation uses asuint64()/asdouble() from math_config.h and arithmetic operations to adjust the exponent, which generates better code than bit masking on ARM and RISC-V architectures. For subnormals, stdc_leading_zeros provides faster normalization than the traditional multiply approach. The zero/infinity/NaN check is simplified to (int64_t)(ix << 1) <= 0, which is more efficient than separate comparisons. Benchmark results on Intel Core i9-13900H (13th Gen): Baseline: 6.778 ns/op Optimized: 4.007 ns/op Speedup: 1.69x (40.9% faster) Zero: 3.580 ns/op (fast path) Denormal: 6.096 ns/op (slower, rare case) Signed-off-by: Osama Abdelkader Reviewed-by: Adhemerval Zanella --- sysdeps/ieee754/dbl-64/s_frexp.c | 55 +++++++++++++------------------- 1 file changed, 22 insertions(+), 33 deletions(-) diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c index 9c45819d4d..7b6818fb78 100644 --- a/sysdeps/ieee754/dbl-64/s_frexp.c +++ b/sysdeps/ieee754/dbl-64/s_frexp.c @@ -18,48 +18,37 @@ #include #include #include +#include +#include "math_config.h" #include -/* - * for non-zero, finite x - * x = frexp(arg,&exp); - * return a double fp quantity x such that 0.5 <= |x| <1.0 - * and the corresponding binary exponent "exp". That is - * arg = x*2^exp. - * If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg - * with *exp=0. - */ - - double __frexp (double x, int *eptr) { - int64_t ix; - EXTRACT_WORDS64 (ix, x); - int32_t ex = 0x7ff & (ix >> 52); - int e = 0; + uint64_t ix = asuint64 (x); + uint32_t ex = (ix >> MANTISSA_WIDTH) & 0x7ff; - 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; - } + int e = ex - EXPONENT_BIAS + 1; + *eptr = e; + return asdouble (ix - ((uint64_t) e << MANTISSA_WIDTH)); + } - ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000); - INSERT_WORDS64 (x, ix); + /* Handle zero, infinity, and NaN. */ + if (__glibc_likely ((int64_t) (ix << 1) <= 0)) + { + *eptr = 0; + return x + x; } - else - /* Quiet signaling NaNs. */ - x += x; - *eptr = e; - return x; + /* Subnormal. */ + uint64_t sign = ix & SIGN_MASK; + int lz = stdc_leading_zeros (ix << (64 - MANTISSA_WIDTH - 1)); + ix <<= lz; + *eptr = -(EXPONENT_BIAS - 2) - lz; + return asdouble ((ix & MANTISSA_MASK) | sign + | (((uint64_t) (EXPONENT_BIAS - 1)) << MANTISSA_WIDTH)); } libm_alias_double (__frexp, frexp) -- 2.47.3