[PATCH v2 3/5] math: Optimize frexp (binary64) with fast path for normal numbers

[Adhemerval Zanella Netto]

On 23/10/25 12:06, Osama Abdelkader wrote:
> Add fast path optimization for frexp using a single unsigned comparison
> to identify normal floating-point numbers and return immediately via
> arithmetic on the bit representation.
> 
> The implementation uses asuint64()/asdouble() 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 (int64_t)(ix << 1) <= 0, which
> is more efficient than separate comparisons.
> 
> Benchmark results on Intel Core i9-13900H (13th Gen):
>   Baseline:     6.778 ns/op
>   Optimized:    4.007 ns/op
>   Speedup:      1.69x (40.9% faster)
>   Zero:         3.580 ns/op (fast path)
>   Denormal:     6.096 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/dbl-64/s_frexp.c | 55 +++++++++++++-------------------
>  1 file changed, 22 insertions(+), 33 deletions(-)
> 
> diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c
> index 9c45819d4d..7b6818fb78 100644
> --- a/sysdeps/ieee754/dbl-64/s_frexp.c
> +++ b/sysdeps/ieee754/dbl-64/s_frexp.c
> @@ -18,48 +18,37 @@
>  #include <inttypes.h>
>  #include <math.h>
>  #include <math_private.h>
> +#include <stdbit.h>
> +#include "math_config.h"
>  #include <libm-alias-double.h>
>  
> -/*
> - * for non-zero, finite x
> - *	x = frexp(arg,&exp);
> - * return a 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 frexp(arg,&exp) returns arg
> - * with *exp=0.
> - */
> -
> -
>  double
>  __frexp (double x, int *eptr)
>  {
> -  int64_t ix;
> -  EXTRACT_WORDS64 (ix, x);
> -  int32_t ex = 0x7ff & (ix >> 52);
> -  int e = 0;
> +  uint64_t ix = asuint64 (x);
> +  uint32_t ex = (ix >> MANTISSA_WIDTH) & 0x7ff;
>  
> -  if (__glibc_likely (ex != 0x7ff && x != 0.0))
> +  /* Fast path for normal numbers.  */
> +  if (__glibc_likely ((ex - 1) < 0x7fe))
>      {
> -      /* Not zero and finite.  */
> -      e = ex - 1022;
> -      if (__glibc_unlikely (ex == 0))
> -	{
> -	  /* Subnormal.  */
> -	  x *= 0x1p54;
> -	  EXTRACT_WORDS64 (ix, x);
> -	  ex = 0x7ff & (ix >> 52);
> -	  e = ex - 1022 - 54;
> -	}
> +      int e = ex - EXPONENT_BIAS + 1;
> +      *eptr = e;
> +      return asdouble (ix - ((uint64_t) e << MANTISSA_WIDTH));
> +    }
>  
> -      ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
> -      INSERT_WORDS64 (x, ix);
> +  /* Handle zero, infinity, and NaN.  */
> +  if (__glibc_likely ((int64_t) (ix << 1) <= 0))
> +    {
> +      *eptr = 0;
> +      return x + x;
>      }
> -  else
> -    /* Quiet signaling NaNs.  */
> -    x += x;
>  
> -  *eptr = e;
> -  return x;
> +  /* Subnormal.  */
> +  uint64_t sign = ix & SIGN_MASK;
> +  int lz = stdc_leading_zeros (ix << (64 - MANTISSA_WIDTH - 1));
> +  ix <<= lz;
> +  *eptr = -(EXPONENT_BIAS - 2) - lz;
> +  return asdouble ((ix & MANTISSA_MASK) | sign
> +		   | (((uint64_t) (EXPONENT_BIAS - 1)) << MANTISSA_WIDTH));
>  }
>  libm_alias_double (__frexp, frexp)