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

[Osama Abdelkader]

Add fast path optimization for frexpl (80-bit x87 extended 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 (se - ex ) to
adjust the exponent directly, which is simpler than bit masking. For subnormals,
the traditional multiply-based normalization is retained as it handles the
split word format more reliably.

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

Benchmark results on Intel Core i9-13900H (13th Gen):
  Baseline:     25.543 ns/op
  Optimized:    25.531 ns/op
  Speedup:      1.00x (neutral)
  Zero:         17.774 ns/op
  Denormal:     23.900 ns/op

Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>
---
v3:
- Address Wilco Dijkstra's review comments.
---
 sysdeps/ieee754/ldbl-96/s_frexpl.c | 98 +++++++++++++++---------------
 1 file changed, 49 insertions(+), 49 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c
index c610704cca..fb4a4d743b 100644
--- a/sysdeps/ieee754/ldbl-96/s_frexpl.c
+++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c
@@ -1,60 +1,60 @@
-/* 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 intel96 (x87 80-bit).
+   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 <float.h>
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-ldouble.h>
 
-static const long double
-#if LDBL_MANT_DIG == 64
-two65 =  3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */
-#else
-# error "Cannot handle this MANT_DIG"
-#endif
+#define EXPONENT_BIAS 16383
 
+static const long double two65 = 0x1p65L;
 
-long double __frexpl(long double x, int *eptr)
+long double
+__frexpl (long double x, int *eptr)
 {
-	uint32_t se, hx, ix, lx;
-	GET_LDOUBLE_WORDS(se,hx,lx,x);
-	ix = 0x7fff&se;
-	*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;
-	SET_LDOUBLE_EXP(x,se);
-	return x;
+  uint32_t se, hx, lx;
+  GET_LDOUBLE_WORDS (se, hx, lx, x);
+  uint32_t ex = se & 0x7fff;
+
+  /* Fast path for normal numbers.  */
+  if (__glibc_likely ((ex - 1) < 0x7ffe))
+    {
+      ex -= EXPONENT_BIAS - 1;
+      *eptr = ex;
+      SET_LDOUBLE_EXP (x, se - ex);
+      return x;
+    }
+
+  /* Handle zero, infinity, and NaN.  */
+  if (__glibc_likely (ex == 0x7fff || ((ex | hx | lx) == 0)))
+    {
+      *eptr = 0;
+      return x + x;
+    }
+
+  /* Subnormal.  */
+  x *= two65;
+  GET_LDOUBLE_EXP (se, x);
+
+  ex = (se & 0x7fff) - EXPONENT_BIAS + 1;
+  *eptr = ex - 65;
+  SET_LDOUBLE_EXP (x, se - ex);
+  return x;
 }
 libm_alias_ldouble (__frexp, frexp)
-- 
2.43.0