math: Optimize frexpf (binary32) with fast path for normal numbers

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>
Reviewed-by: Adhemerval Zanella  <adhemerval.zanella@linaro.org>
This commit is contained in:
Osama Abdelkader 2025-10-23 18:06:28 +03:00 committed by Adhemerval Zanella
parent ff041e8f8e
commit 4d2582150e
1 changed files with 44 additions and 34 deletions

View File

@ -1,44 +1,54 @@
/* s_frexpf.c -- float version of s_frexp.c.
*/
/* Optimized frexp implementation for binary32.
Copyright (C) 2025 Free Software Foundation, Inc.
This file is part of the GNU C Library.
/*
* ====================================================
* 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.
* ====================================================
*/
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.
#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
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)