[PATCH 2/5] math: Optimize frexpf (binary32) with fast path for normal numbers

[Osama Abdelkader]

Add fast path optimization for frexpf using a single unsigned comparison
to identify normal floating-point numbers and return immediately via bit
manipulation.

The implementation uses asuint()/asfloat() from math_config.h and symbolic
constants (MANTISSA_MASK, SIGN_MASK, EXPONENT_BIAS, MANTISSA_WIDTH) for
better readability and maintainability.

For normal numbers (exponent 1 to 254), the function extracts the unbiased
exponent and constructs the result in [0.5, 1.0) using direct bit operations,
avoiding branches and floating-point operations.

Benchmark results on Intel Core i9-13900H (13th Gen):
  Baseline:  5.858 ns/op
  Optimized: 5.878 ns/op (average of 3 runs)
  Speedup:   Neutral (within measurement noise)

Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>
---
 sysdeps/ieee754/flt-32/s_frexpf.c | 37 +++++++++++++++++++------------
 1 file changed, 23 insertions(+), 14 deletions(-)

diff --git a/sysdeps/ieee754/flt-32/s_frexpf.c b/sysdeps/ieee754/flt-32/s_frexpf.c
index 59fef66a6b..da08c9367c 100644
--- a/sysdeps/ieee754/flt-32/s_frexpf.c
+++ b/sysdeps/ieee754/flt-32/s_frexpf.c
@@ -18,27 +18,36 @@ 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
 two25 =  3.3554432000e+07; /* 0x4c000000 */
 
+/* Exponent bits for result in [0.5, 1.0): 0x7e = -1 + EXPONENT_BIAS.  */
+#define FREXPF_EXP \
+  (((uint32_t) (EXPONENT_BIAS - 1)) << MANTISSA_WIDTH)
+
 float __frexpf(float x, int *eptr)
 {
-	int32_t hx,ix;
-	GET_FLOAT_WORD(hx,x);
-	ix = 0x7fffffff&hx;
+	uint32_t hx = asuint (x);
+	uint32_t ex = (hx >> MANTISSA_WIDTH) & 0xff;
+
+	/* Fast path for normal numbers.  */
+	if (__glibc_likely ((ex - 1) < 0xfe))
+	  {
+	    *eptr = ex - EXPONENT_BIAS + 1;
+	    return asfloat ((hx & (MANTISSA_MASK | SIGN_MASK)) | FREXPF_EXP);
+	  }
+
+	/* Handle special cases: zero, subnormal, infinity, NaN.  */
+	uint32_t ix = hx & 0x7fffffff;
 	*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;
-	SET_FLOAT_WORD(x,hx);
-	return x;
+	if (ix >= 0x7f800000 || (ix == 0)) return x + x;	/* 0,inf,nan */
+	/* Subnormal */
+	x *= two25;
+	hx = asuint (x);
+	*eptr = ((hx & 0x7fffffff) >> MANTISSA_WIDTH) - EXPONENT_BIAS - 25 + 1;
+	return asfloat ((hx & (MANTISSA_MASK | SIGN_MASK)) | FREXPF_EXP);
 }
 libm_alias_float (__frexp, frexp)
-- 
2.43.0