[PATCH v2 5/5] math: Optimize frexpl (binary128) with fast path for normal numbers

[Adhemerval Zanella Netto]

On 23/10/25 12:06, Osama Abdelkader wrote:
> Add fast path optimization for frexpl (128-bit IEEE quad 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 (hx - ((ex - FREXPL_EXP_VALUE) << 48))
> to adjust the exponent in place, which is simpler and more efficient than
> bit masking. For subnormals, the traditional multiply-based normalization
> is retained for reliability with the split 64-bit word format.
> 
> The zero/infinity/NaN check groups these special cases together for better
> branch prediction.
> 
> This optimization provides the same algorithmic improvements as the other
> frexp variants while maintaining correctness for all edge cases.
> 
> Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>

LGTM, thanks.  On aarch64 N1 I see:

master:

   "workload-random-m20-p20": {
    "duration": 1.04133e+09,
    "iterations": 7.2e+07,
    "reciprocal-throughput": 5.31921,
    "latency": 23.6066,
    "max-throughput": 1.87998e+08,
    "min-throughput": 4.2361e+07
   }
  }

patched:

   "workload-random-m20-p20": {
    "duration": 1.04022e+09,
    "iterations": 7.6e+07,
    "reciprocal-throughput": 4.72723,
    "latency": 22.6469,
    "max-throughput": 2.1154e+08,
    "min-throughput": 4.41561e+07
   }
  }

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


> ---
>  sysdeps/ieee754/ldbl-128/s_frexpl.c | 100 +++++++++++++++-------------
>  1 file changed, 55 insertions(+), 45 deletions(-)
> 
> diff --git a/sysdeps/ieee754/ldbl-128/s_frexpl.c b/sysdeps/ieee754/ldbl-128/s_frexpl.c
> index e4db093a05..0e56e86878 100644
> --- a/sysdeps/ieee754/ldbl-128/s_frexpl.c
> +++ b/sysdeps/ieee754/ldbl-128/s_frexpl.c
> @@ -1,54 +1,64 @@
> -/* 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 binary128 (IEEE quad precision).
> +   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 <libm-alias-ldouble.h>
>  
> -static const _Float128
> -two114 = L(2.0769187434139310514121985316880384E+34); /* 0x4071000000000000, 0 */
> +/* Exponent value for result in [0.5, 1.0): 0x3ffe = -1 + 16383.  */
> +#define FREXPL_EXP_VALUE 0x3ffe
> +#define EXPONENT_BIAS 16383
> +#define MANTISSA_MASK UINT64_C (0x0000ffffffffffff)
> +#define SIGN_MASK UINT64_C (0x8000000000000000)
> +
> +static const _Float128 two114 = 0x1p114L;
>  
> -_Float128 __frexpl(_Float128 x, int *eptr)
> +_Float128
> +__frexpl (_Float128 x, int *eptr)
>  {
> -	uint64_t hx, lx, ix;
> -	GET_LDOUBLE_WORDS64(hx,lx,x);
> -	ix = 0x7fffffffffffffffULL&hx;
> -	*eptr = 0;
> -	if(ix>=0x7fff000000000000ULL||((ix|lx)==0)) return x + x;/* 0,inf,nan */
> -	if (ix<0x0001000000000000ULL) {		/* subnormal */
> -	    x *= two114;
> -	    GET_LDOUBLE_MSW64(hx,x);
> -	    ix = hx&0x7fffffffffffffffULL;
> -	    *eptr = -114;
> -	}
> -	*eptr += (ix>>48)-16382;
> -	hx = (hx&0x8000ffffffffffffULL) | 0x3ffe000000000000ULL;
> -	SET_LDOUBLE_MSW64(x,hx);
> -	return x;
> +  uint64_t hx, lx;
> +  GET_LDOUBLE_WORDS64 (hx, lx, x);
> +  uint64_t ex = 0x7fff & (hx >> 48);
> +
> +  /* Fast path for normal numbers.  */
> +  if (__glibc_likely ((ex - 1U) < 0x7ffe))
> +    {
> +      *eptr = ex - EXPONENT_BIAS + 1;
> +      hx = hx - ((ex - FREXPL_EXP_VALUE) << 48);
> +      SET_LDOUBLE_MSW64 (x, hx);
> +      return x;
> +    }
> +
> +  /* Handle zero, infinity, and NaN.  */
> +  uint64_t ix = hx & 0x7fffffffffffffffULL;
> +  if (__glibc_likely (ix >= 0x7fff000000000000ULL || ((ix | lx) == 0)))
> +    {
> +      *eptr = 0;
> +      return x + x;
> +    }
> +
> +  /* Subnormal.  */
> +  x *= two114;
> +  GET_LDOUBLE_MSW64 (hx, x);
> +  ex = 0x7fff & (hx >> 48);
> +  *eptr = ex - EXPONENT_BIAS - 114 + 1;
> +  hx = hx - ((ex - FREXPL_EXP_VALUE) << 48);
> +  SET_LDOUBLE_MSW64 (x, hx);
> +  return x;
>  }
>  libm_alias_ldouble (__frexp, frexp)