[PATCH v2 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
arithmetic on the bit representation.

The implementation uses asuint()/asfloat() from math_config.h and arithmetic
operations to adjust the exponent, which generates better code than bit masking
on ARM and RISC-V architectures. For subnormals, stdc_leading_zeros provides
faster normalization than the traditional multiply approach.

The zero/infinity/NaN check is simplified to (int32_t)(hx << 1) <= 0, which
is more efficient than separate comparisons.

Benchmark results on Intel Core i9-13900H (13th Gen):
  Baseline:     5.858 ns/op
  Optimized:    4.003 ns/op
  Speedup:      1.46x (31.7% faster)
  Zero:         3.580 ns/op (fast path)
  Denormal:     5.597 ns/op (slower, rare case)

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

diff --git a/sysdeps/ieee754/flt-32/s_frexpf.c b/sysdeps/ieee754/flt-32/s_frexpf.c
index 59fef66a6b..f424e52199 100644
--- a/sysdeps/ieee754/flt-32/s_frexpf.c
+++ b/sysdeps/ieee754/flt-32/s_frexpf.c
@@ -1,44 +1,54 @@
-/* s_frexpf.c -- float 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: s_frexpf.c,v 1.5 1995/05/10 20:47:26 jtc Exp $";
-#endif
+/* Optimized frexp implementation for binary32.
+   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 <stdbit.h>
+#include "math_config.h"
 #include <libm-alias-float.h>
 
-static const float
-two25 =  3.3554432000e+07; /* 0x4c000000 */
-
-float __frexpf(float x, int *eptr)
+float
+__frexpf (float x, int *eptr)
 {
-	int32_t hx,ix;
-	GET_FLOAT_WORD(hx,x);
-	ix = 0x7fffffff&hx;
-	*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;
+  uint32_t hx = asuint (x);
+  uint32_t ex = (hx >> MANTISSA_WIDTH) & 0xff;
+
+  /* Fast path for normal numbers.  */
+  if (__glibc_likely ((ex - 1) < 0xfe))
+    {
+      int e = ex - EXPONENT_BIAS + 1;
+      *eptr = e;
+      return asfloat (hx - (e << MANTISSA_WIDTH));
+    }
+
+  /* Handle zero, infinity, and NaN.  */
+  if (__glibc_likely ((int32_t) (hx << 1) <= 0))
+    {
+      *eptr = 0;
+      return x + x;
+    }
+
+  /* Subnormal.  */
+  uint32_t sign = hx & SIGN_MASK;
+  int lz = stdc_leading_zeros (hx << (32 - MANTISSA_WIDTH - 1));
+  hx <<= lz;
+  *eptr = -(EXPONENT_BIAS - 2) - lz;
+  return asfloat ((hx & MANTISSA_MASK) | sign
+		  | (((uint32_t) (EXPONENT_BIAS - 1)) << MANTISSA_WIDTH));
 }
 libm_alias_float (__frexp, frexp)
-- 
2.43.0