math: Optimize frexp (binary64) with fast path for normal numbers

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

View File

@ -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;
}
ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
INSERT_WORDS64 (x, ix);
int e = ex - EXPONENT_BIAS + 1;
*eptr = e;
return asdouble (ix - ((uint64_t) e << MANTISSA_WIDTH));
}
else
/* Quiet signaling NaNs. */
x += x;
*eptr = e;
return x;
/* Handle zero, infinity, and NaN. */
if (__glibc_likely ((int64_t) (ix << 1) <= 0))
{
*eptr = 0;
return x + 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)