[PATCH v2 2/5] math: Optimize frexpf (binary32) with fast path for normal numbers

[Adhemerval Zanella Netto]

On 23/10/25 12:06, Osama Abdelkader wrote:
> 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>


LGTM, thanks.  I am seeing similar numbers on aarch64 N1.

Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>

> ---
>  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)