From 4f18501498e8fe7eca8d4adaeabb108ae62972cd Mon Sep 17 00:00:00 2001 From: Osama Abdelkader Date: Mon, 10 Nov 2025 20:17:07 +0000 Subject: [PATCH] math: Optimize frexpl (intel96) with fast path for normal numbers Add fast path optimization for frexpl (80-bit x87 extended 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 (se - ex ) to adjust the exponent directly, which is simpler than bit masking. For subnormals, the traditional multiply-based normalization is retained as it handles the split word format more reliably. The zero/infinity/NaN check groups these special cases together for better branch prediction. Benchmark results on Intel Core i9-13900H (13th Gen): Baseline: 25.543 ns/op Optimized: 25.531 ns/op Speedup: 1.00x (neutral) Zero: 17.774 ns/op Denormal: 23.900 ns/op Signed-off-by: Osama Abdelkader Reviewed-by: Wilco Dijkstra --- sysdeps/ieee754/ldbl-96/s_frexpl.c | 98 +++++++++++++++--------------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c index c610704cca9..c92dcbedcc8 100644 --- a/sysdeps/ieee754/ldbl-96/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c @@ -1,60 +1,60 @@ -/* 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 intel96 (x87 80-bit). + 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 + . */ #include #include #include #include -static const long double -#if LDBL_MANT_DIG == 64 -two65 = 3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */ -#else -# error "Cannot handle this MANT_DIG" -#endif +#define EXPONENT_BIAS 16383 +static const long double two65 = 0x1p65L; -long double __frexpl(long double x, int *eptr) +long double +__frexpl (long double x, int *eptr) { - uint32_t se, hx, ix, lx; - GET_LDOUBLE_WORDS(se,hx,lx,x); - ix = 0x7fff&se; - *eptr = 0; - if(ix==0x7fff||((ix|hx|lx)==0)) return x + x; /* 0,inf,nan */ - if (ix==0x0000) { /* subnormal */ - x *= two65; - GET_LDOUBLE_EXP(se,x); - ix = se&0x7fff; - *eptr = -65; - } - *eptr += ix-16382; - se = (se & 0x8000) | 0x3ffe; - SET_LDOUBLE_EXP(x,se); - return x; + uint32_t se, hx, lx; + GET_LDOUBLE_WORDS (se, hx, lx, x); + uint32_t ex = se & 0x7fff; + + /* Fast path for normal numbers. */ + if (__glibc_likely ((ex - 1) < 0x7ffe)) + { + ex -= EXPONENT_BIAS - 1; + *eptr = ex; + SET_LDOUBLE_EXP (x, se - ex); + return x; + } + + /* Handle zero, infinity, and NaN. */ + if (__glibc_likely (ex == 0x7fff || (ex | hx | lx) == 0)) + { + *eptr = 0; + return x + x; + } + + /* Subnormal. */ + x *= two65; + GET_LDOUBLE_EXP (se, x); + + ex = (se & 0x7fff) - EXPONENT_BIAS + 1; + *eptr = ex - 65; + SET_LDOUBLE_EXP (x, se - ex); + return x; } libm_alias_ldouble (__frexp, frexp) -- 2.47.3