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

[Osama Abdelkader]

Add fast path optimization for frexp using a single unsigned comparison
to identify normal floating-point numbers and return immediately via bit
manipulation.

The implementation uses asuint64()/asdouble() from math_config.h and symbolic
constants (MANTISSA_MASK, SIGN_MASK, EXPONENT_BIAS, MANTISSA_WIDTH) for
better readability and maintainability.

For normal numbers (exponent 1 to 2046), the function extracts the unbiased
exponent and constructs the result in [0.5, 1.0) using direct bit operations,
avoiding branches and floating-point operations.

Benchmark results on Intel Core i9-13900H (13th Gen):
  Baseline:  6.778 ns/op
  Optimized: 6.210 ns/op (average of 3 runs)
  Speedup:   1.09x (8.4% faster)

Signed-off-by: Osama Abdelkader <osama.abdelkader@gmail.com>
---
 sysdeps/ieee754/dbl-64/s_frexp.c | 39 ++++++++++++++++++--------------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c
index 9c45819d4d..08a5c5575e 100644
--- a/sysdeps/ieee754/dbl-64/s_frexp.c
+++ b/sysdeps/ieee754/dbl-64/s_frexp.c
@@ -18,6 +18,7 @@
 #include <inttypes.h>
 #include <math.h>
 #include <math_private.h>
+#include "math_config.h"
 #include <libm-alias-double.h>
 
 /*
@@ -30,30 +31,34 @@
  * with *exp=0.
  */
 
+/* Exponent bits for result in [0.5, 1.0): 0x3fe = -1 + EXPONENT_BIAS.  */
+#define FREXP_EXP \
+  (((uint64_t) (EXPONENT_BIAS - 1)) << MANTISSA_WIDTH)
 
 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;
-	}
+      *eptr = ex - EXPONENT_BIAS + 1;
+      return asdouble ((ix & (MANTISSA_MASK | SIGN_MASK)) | FREXP_EXP);
+    }
 
-      ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000);
-      INSERT_WORDS64 (x, ix);
+  /* Handle special cases: zero, subnormal, infinity, NaN.  */
+  int e = 0;
+  if (__glibc_likely (ex != 0x7ff && x != 0.0))
+    {
+      /* Subnormal.  */
+      x *= 0x1p54;
+      ix = asuint64 (x);
+      ex = (ix >> MANTISSA_WIDTH) & 0x7ff;
+      e = ex - EXPONENT_BIAS - 54 + 1;
+      *eptr = e;
+      return asdouble ((ix & (MANTISSA_MASK | SIGN_MASK)) | FREXP_EXP);
     }
   else
     /* Quiet signaling NaNs.  */
-- 
2.43.0