math: Simplify and optimize modf implementation

Refactor the generic implementation to use math_config.h definitions,
and add an alternative one if the ABI supports truncf instructions
(gated through math-use-builtins-trunc.h).

The generic implementation generates similar code on x86_64, while
the optimization one for aarch64 (where truncf is supported as a
builtin by through frintz), the improvements are:

reciprocal-throughput           master    patch    difference
workload-0_1                    3.0595   3.0698        -0.34%
workload-1_maxint               5.1747   3.0542        40.98%
workload-maxint_maxfloat        3.4391   3.0349        11.75%
workload-integral               3.2732   3.0293         7.45%

latency                         master    patch    difference
workload-0_1                    3.5267   4.7107       -33.57%
workload-1_maxint               6.9074   4.7282        31.55%
workload-maxint_maxfloat        3.7210   4.7506       -27.67%
workload-integral               3.8634   4.8137       -24.60%

Checked on aarch64-linux-gnu and x86_64-linux-gnu.
Reviewed-by: Wilco Dijkstra  <Wilco.Dijkstra@arm.com>

diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h
index 5766f6a1db91911d2249a73f4f8d075f28f068f6..d9288c4fb391aec3bb70409439ab268f38519fe6 100644 (file)
--- a/sysdeps/ieee754/dbl-64/math_config.h
+++ b/sysdeps/ieee754/dbl-64/math_config.h
@@ -109,6 +109,7 @@ issignaling_inline (double x)
 #define BIT_WIDTH       64
 #define MANTISSA_WIDTH  52
 #define EXPONENT_WIDTH  11
+#define EXPONENT_BIAS   1023
 #define MANTISSA_MASK   UINT64_C(0x000fffffffffffff)
 #define EXPONENT_MASK   UINT64_C(0x7ff0000000000000)
 #define EXP_MANT_MASK   UINT64_C(0x7fffffffffffffff)
@@ -121,12 +122,24 @@ is_nan (uint64_t x)
   return (x & EXP_MANT_MASK) > EXPONENT_MASK;
 }
 
+static inline bool
+is_inf (uint64_t x)
+{
+  return (x << 1) == (EXPONENT_MASK << 1);
+}
+
 static inline uint64_t
 get_mantissa (uint64_t x)
 {
   return x & MANTISSA_MASK;
 }
 
+static inline int
+get_exponent (uint64_t x)
+{
+  return (int)((x >> MANTISSA_WIDTH & 0x7ff) - EXPONENT_BIAS);
+}
+
 /* Convert integer number X, unbiased exponent EP, and sign S to double:
 
    result = X * 2^(EP+1 - exponent_bias)

diff --git a/sysdeps/ieee754/dbl-64/s_modf.c b/sysdeps/ieee754/dbl-64/s_modf.c
index 0de2084caf11dc8eb522c6805d7982ec03dc5a9e..90cd8e8c3ef8db2c925da074f4e0b48e96c15afa 100644 (file)
--- a/sysdeps/ieee754/dbl-64/s_modf.c
+++ b/sysdeps/ieee754/dbl-64/s_modf.c
@@ -1,63 +1,68 @@
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
+/* Extract signed integral and fractional values.
+   Copyright (C) 1993-2025 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
 
-/*
- * modf(double x, double *iptr)
- * return fraction part of x, and return x's integral part in *iptr.
- * Method:
- *     Bit twiddling.
- *
- * Exception:
- *     No exception.
- */
+   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-double.h>
-#include <stdint.h>
-
-static const double one = 1.0;
+#include "math_config.h"
+#include <math-use-builtins-trunc.h>
 
 double
-__modf(double x, double *iptr)
+__modf (double x, double *iptr)
 {
-       int64_t i0;
-       int32_t j0;
-       EXTRACT_WORDS64(i0,x);
-       j0 = ((i0>>52)&0x7ff)-0x3ff;    /* exponent of x */
-       if(j0<52) {                     /* integer part in x */
-           if(j0<0) {                  /* |x|<1 */
-               /* *iptr = +-0 */
-               INSERT_WORDS64(*iptr,i0&UINT64_C(0x8000000000000000));
-               return x;
-           } else {
-               uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
-               if((i0&i)==0) {         /* x is integral */
-                   *iptr = x;
-                   /* return +-0 */
-                   INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));
-                   return x;
-               } else {
-                   INSERT_WORDS64(*iptr,i0&(~i));
-                   return x - *iptr;
-               }
-           }
-       } else { /* no fraction part */
-           *iptr = x*one;
-           /* We must handle NaNs separately.  */
-           if (j0 == 0x400 && (i0 & UINT64_C(0xfffffffffffff)))
-             return x*one;
-           INSERT_WORDS64(x,i0&UINT64_C(0x8000000000000000));  /* return +-0 */
-           return x;
+  uint64_t t = asuint64 (x);
+#if USE_TRUNC_BUILTIN
+  if (is_inf (t))
+    {
+      *iptr = x;
+      return copysign (0.0, x);
+    }
+  *iptr = trunc (x);
+  return copysign (x - *iptr, x);
+#else
+  int e = get_exponent (t);
+  /* No fraction part.  */
+  if (e < MANTISSA_WIDTH)
+    {
+      if (e < 0)
+       {
+         /* |x|<1 -> *iptr = +-0 */
+         *iptr = asdouble (t & SIGN_MASK);
+         return x;
+       }
+
+      uint64_t i = MANTISSA_MASK >> e;
+      if ((t & i) == 0)
+       {
+         /* x in integral, return +-0  */
+         *iptr = x;
+         return asdouble (t & SIGN_MASK);
        }
+
+      *iptr = asdouble (t & ~i);
+      return x - *iptr;
+    }
+
+  /* Set invalid operation for sNaN.  */
+  *iptr = x * 1.0;
+  if ((e == 0x400) && (t & MANTISSA_MASK))
+    return *iptr;
+  return asdouble (t & SIGN_MASK);
+#endif
 }
 #ifndef __modf
 libm_alias_double (__modf, modf)