[PATCH 4/5] math: Optimize frexpl (intel96) with fast path for normal numbers

[Osama Abdelkader]

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

The implementation uses symbolic constants (FREXPL_EXP, EXPONENT_BIAS,
SIGN_MASK) for better readability and maintainability.

For normal numbers (exponent 1 to 32766), the function extracts the unbiased
exponent and sets the result exponent to -1 + bias (0x3ffe) to place the
value in [0.5, 1.0), avoiding branches and extra operations.

Benchmark results on Intel Core i9-13900H (13th Gen):
  Baseline:  25.543 ns/op
  Optimized: 24.402 ns/op (average of 3 runs)
  Speedup:   1.05x (4.5% faster)

Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>
---
 sysdeps/ieee754/ldbl-96/s_frexpl.c | 35 ++++++++++++++++++++----------
 1 file changed, 24 insertions(+), 11 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c
index c610704cca..2ce3d7d9e4 100644
--- a/sysdeps/ieee754/ldbl-96/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c
@@ -38,22 +38,35 @@ two65 =  3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */
 # error "Cannot handle this MANT_DIG"
 #endif
 
+/* Biased exponent for result in [0.5, 1.0): 0x3ffe = -1 + 16383.  */
+#define FREXPL_EXP 0x3ffe
+#define EXPONENT_BIAS 16383
+#define SIGN_MASK 0x8000
 
 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 - EXPONENT_BIAS + 1;
+	    se = (se & SIGN_MASK) | FREXPL_EXP;
+	    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;
-	se = (se & 0x8000) | 0x3ffe;
+	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 - EXPONENT_BIAS - 65 + 1;
+	se = (se & SIGN_MASK) | FREXPL_EXP;
 	SET_LDOUBLE_EXP(x,se);
 	return x;
 }
-- 
2.43.0