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

[Osama Abdelkader]

Add fast path optimization for frexpl (128-bit IEEE quad precision) 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,
MANTISSA_MASK, SIGN_MASK) for better readability and maintainability.

For normal numbers (exponent 1 to 32766), the function extracts the unbiased
exponent and constructs the result in [0.5, 1.0) by setting the biased
exponent to -1 + 16383 (0x3ffe), avoiding branches and extra operations.

This optimization provides the same fast path benefits as the other frexp
variants, with cleaner code using symbolic constants instead of magic
numbers.

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

diff --git a/sysdeps/ieee754/ldbl-128/s_frexpl.c b/sysdeps/ieee754/ldbl-128/s_frexpl.c
index e4db093a05..9f18c2fee0 100644
--- a/sysdeps/ieee754/ldbl-128/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-128/s_frexpl.c
@@ -33,21 +33,37 @@ static char rcsid[] = "$NetBSD: $";
 static const _Float128
 two114 = L(2.0769187434139310514121985316880384E+34); /* 0x4071000000000000, 0 */
 
+/* Biased exponent for result in [0.5, 1.0): 0x3ffe = -1 + 16383.  */
+#define FREXPL_EXP UINT64_C(0x3ffe000000000000)
+#define EXPONENT_BIAS 16383
+#define MANTISSA_MASK UINT64_C(0x0000ffffffffffff)
+#define SIGN_MASK UINT64_C(0x8000000000000000)
+
 _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 - EXPONENT_BIAS + 1;
+	    hx = (hx & (SIGN_MASK | MANTISSA_MASK)) | FREXPL_EXP;
+	    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);
+	ex = 0x7fff & (hx >> 48);
+	*eptr = ex - EXPONENT_BIAS - 114 + 1;
+	hx = (hx & (SIGN_MASK | MANTISSA_MASK)) | FREXPL_EXP;
 	SET_LDOUBLE_MSW64(x,hx);
 	return x;
 }
-- 
2.43.0