[PATCH v2 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 arithmetic on the exponent field.

The implementation uses arithmetic operations (hx - ((ex - FREXPL_EXP_VALUE) << 48))
to adjust the exponent in place, which is simpler and more efficient than
bit masking. For subnormals, the traditional multiply-based normalization
is retained for reliability with the split 64-bit word format.

The zero/infinity/NaN check groups these special cases together for better
branch prediction.

This optimization provides the same algorithmic improvements as the other
frexp variants while maintaining correctness for all edge cases.

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

diff --git a/sysdeps/ieee754/ldbl-128/s_frexpl.c b/sysdeps/ieee754/ldbl-128/s_frexpl.c
index e4db093a05..0e56e86878 100644
--- a/sysdeps/ieee754/ldbl-128/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-128/s_frexpl.c
@@ -1,54 +1,64 @@
-/* s_frexpl.c -- long double version of s_frexp.c.
- */
-
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: $";
-#endif
-
-/*
- * for non-zero x
- *	x = frexpl(arg,&exp);
- * return a long 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 frexpl(arg,&exp) returns arg
- * with *exp=0.
- */
+/* Optimized frexp implementation for binary128 (IEEE quad precision).
+   Copyright (C) 2025 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
 
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-ldouble.h>
 
-static const _Float128
-two114 = L(2.0769187434139310514121985316880384E+34); /* 0x4071000000000000, 0 */
+/* Exponent value for result in [0.5, 1.0): 0x3ffe = -1 + 16383.  */
+#define FREXPL_EXP_VALUE 0x3ffe
+#define EXPONENT_BIAS 16383
+#define MANTISSA_MASK UINT64_C (0x0000ffffffffffff)
+#define SIGN_MASK UINT64_C (0x8000000000000000)
+
+static const _Float128 two114 = 0x1p114L;
 
-_Float128 __frexpl(_Float128 x, int *eptr)
+_Float128
+__frexpl (_Float128 x, int *eptr)
 {
-	uint64_t hx, lx, ix;
-	GET_LDOUBLE_WORDS64(hx,lx,x);
-	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;
-	SET_LDOUBLE_MSW64(x,hx);
-	return x;
+  uint64_t hx, lx;
+  GET_LDOUBLE_WORDS64 (hx, lx, x);
+  uint64_t ex = 0x7fff & (hx >> 48);
+
+  /* Fast path for normal numbers.  */
+  if (__glibc_likely ((ex - 1U) < 0x7ffe))
+    {
+      *eptr = ex - EXPONENT_BIAS + 1;
+      hx = hx - ((ex - FREXPL_EXP_VALUE) << 48);
+      SET_LDOUBLE_MSW64 (x, hx);
+      return x;
+    }
+
+  /* Handle zero, infinity, and NaN.  */
+  uint64_t ix = hx & 0x7fffffffffffffffULL;
+  if (__glibc_likely (ix >= 0x7fff000000000000ULL || ((ix | lx) == 0)))
+    {
+      *eptr = 0;
+      return x + x;
+    }
+
+  /* Subnormal.  */
+  x *= two114;
+  GET_LDOUBLE_MSW64 (hx, x);
+  ex = 0x7fff & (hx >> 48);
+  *eptr = ex - EXPONENT_BIAS - 114 + 1;
+  hx = hx - ((ex - FREXPL_EXP_VALUE) << 48);
+  SET_LDOUBLE_MSW64 (x, hx);
+  return x;
 }
 libm_alias_ldouble (__frexp, frexp)
-- 
2.43.0