math: Remove the SVID error handling from remainder

The optimized i386 version is faster than the generic one, and gcc implements it through the builtin. This optimization enables us to migrate the implementation to a C version. The performance on a Zen3 chip is similar to the SVID one. The m68k provided an optimized version through __m81_u(remainderf) (mathimpl.h), and gcc does not implement it through a builtin (different than i386). Performance improves a bit on x86_64 (Zen3, gcc 15.2.1): reciprocal-throughput input master NO-SVID improvement x86_64 subnormals 18.8522 16.2506 13.80% x86_64 normal 421.8260 403.9270 4.24% x86_64 close-exponent 21.0579 18.7642 10.89% i686 subnormals 21.3443 21.4229 -0.37% i686 normal 525.8380 538.807 -2.47% i686 close-exponent 21.6589 21.7983 -0.64% Tested on x86_64-linux-gnu and i686-linux-gnu. Reviewed-by: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
2025-11-15 15:21:18 +03:00 · 2025-10-31 13:08:51 -03:00
parent c4c6c79d70
commit 0e1a1178ee
33 changed files with 118 additions and 37 deletions
--- a/math/Versions
+++ b/math/Versions
@@ -697,6 +697,7 @@ libm {
    j1f;
    jnf;
    log10f;
+    remainder;
    remainderf;
    y0f;
    y1f;
--- a/math/w_remainder_compat.c
+++ b/math/w_remainder_compat.c
@@ -19,23 +19,32 @@
 #include <math_private.h>
 #include <math-svid-compat.h>
 #include <libm-alias-double.h>
+#include <shlib-compat.h>


-#if LIBM_SVID_COMPAT
+#if LIBM_SVID_COMPAT && (SHLIB_COMPAT (libm, GLIBC_2_0, GLIBC_2_43) \
+			 || defined NO_LONG_DOUBLE \
+			 || defined LONG_DOUBLE_COMPAT)
 /* wrapper remainder */
 double
-__remainder (double x, double y)
+__remainder_compat (double x, double y)
 {
  if (((__builtin_expect (y == 0.0, 0) && ! isnan (x))
       || (__builtin_expect (isinf (x), 0) && ! isnan (y)))
      && _LIB_VERSION != _IEEE_)
    return __kernel_standard (x, y, 28); /* remainder domain */

-  return __ieee754_remainder (x, y);
+  return __remainder (x, y);
 }
-libm_alias_double (__remainder, remainder)
-weak_alias (__remainder, drem)
+compat_symbol (libm, __remainder_compat, remainder, GLIBC_2_0);
+weak_alias (__remainder_compat, drem)
 # ifdef NO_LONG_DOUBLE
-weak_alias (__remainder, dreml)
+weak_alias (__remainder_compat, dreml)
+weak_alias (__remainder_compat, remainderl)
+# endif
+# ifdef LONG_DOUBLE_COMPAT
+LONG_DOUBLE_COMPAT_CHOOSE_libm_remainderl (
+  compat_symbol (libm, __remainder_compat, remainderl, \
+		 FIRST_VERSION_libm_remainderl), );
 # endif
 #endif
--- a/sysdeps/i386/fpu/e_remainder.S
+++ b/sysdeps/i386/fpu/e_remainder.S
@@ -1,18 +0,0 @@
-/*
- * Public domain.
- */
-
-#include <machine/asm.h>
-#include <libm-alias-finite.h>
-
-ENTRY(__ieee754_remainder)
-	fldl	12(%esp)
-	fldl	4(%esp)
-1:	fprem1
-	fstsw	%ax
-	sahf
-	jp	1b
-	fstp	%st(1)
-	ret
-END (__ieee754_remainder)
-libm_alias_finite (__ieee754_remainder, __remainder)
--- a/sysdeps/i386/fpu/e_remainder.c
+++ b/sysdeps/i386/fpu/e_remainder.c
@@ -0,0 +1,41 @@
+/* Floating-point remainder function.
+   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 <libm-alias-finite.h>
+#include <libm-alias-double.h>
+#include "math_config.h"
+
+double
+__remainder (double x, double y)
+{
+  uint64_t hx = asuint64 (x);
+  uint64_t hy = asuint64 (y);
+
+  /* fmod(+-Inf,y) or fmod(x,0) */
+  if (__glibc_unlikely ((is_inf (hx) || y == 0.0)
+			&& !is_nan (hy)
+			&& !is_nan (hx)))
+    return __math_invalid (x);
+
+  return __builtin_remainder (x, y);
+}
+strong_alias (__remainder, __ieee754_remainder)
+versioned_symbol (libm, __remainder, remainder, GLIBC_2_43);
+libm_alias_double_other (__remainder, remainder)
+libm_alias_finite (__ieee754_remainder, __remainder)
--- a/sysdeps/ieee754/dbl-64/e_remainder.c
+++ b/sysdeps/ieee754/dbl-64/e_remainder.c
@@ -18,10 +18,12 @@

 #include <math.h>
 #include <libm-alias-finite.h>
+#include <libm-alias-double.h>
+#include <math-svid-compat.h>
 #include "math_config.h"

 double
-__ieee754_remainder (double x, double y)
+__remainder (double x, double y)
 {
  uint64_t hx = asuint64 (x);
  uint64_t hy = asuint64 (y);
@@ -34,12 +36,8 @@ __ieee754_remainder (double x, double y)
  y = fabs (y);
  if (__glibc_likely (hy < UINT64_C (0x7fe0000000000000)))
    {
-      /* |x| not finite, |y| equal 0 is handled by fmod.  */
-      if (__glibc_unlikely (hx >= EXPONENT_MASK))
-	return (x * y) / (x * y);
-
-      x = fabs (__ieee754_fmod (x, y + y));
-      if (x + x > y)
+      x = fabs (__fmod (x, y + y));
+      if (isgreater (x + x, y))
 	{
 	  x -= y;
 	  if (x + x >= y)
@@ -52,9 +50,9 @@ __ieee754_remainder (double x, double y)
    }
  else
    {
-      /* |x| not finite or |y| is NaN or 0 */
-      if ((hx >= EXPONENT_MASK || (hy - 1) >= EXPONENT_MASK))
-	return (x * y) / (x * y);
+      /* |x| not finite or |y| is NaN */
+      if (__glibc_unlikely (hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
+	return __math_invalid (x * y);

      x = fabs (x);
      double y_half = y * 0.5;
@@ -70,4 +68,14 @@ __ieee754_remainder (double x, double y)

  return sx ? -x : x;
 }
-libm_alias_finite (__ieee754_remainder, __remainder)
+libm_alias_finite (__remainder, __remainder)
+#if LIBM_SVID_COMPAT
+versioned_symbol (libm, __remainder, remainder, GLIBC_2_43);
+libm_alias_double_other (__remainder, remainder)
+#else
+libm_alias_double (__remainder, remainder)
+weak_alias (__remainder, drem)
+# ifdef NO_LONG_DOUBLE
+weak_alias (__remainder, dreml)
+# endif
+#endif
--- a/sysdeps/ieee754/dbl-64/w_remainder.c
+++ b/sysdeps/ieee754/dbl-64/w_remainder.c
@@ -0,0 +1 @@
+/* Not needed */
--- a/sysdeps/ieee754/ldbl-opt/w_remainder_compat.c
+++ b/sysdeps/ieee754/ldbl-opt/w_remainder_compat.c
@@ -1,6 +1,6 @@
 #include <math_ldbl_opt.h>
 #include <math/w_remainder_compat.c>
 #if LONG_DOUBLE_COMPAT(libm, GLIBC_2_0)
-strong_alias (__remainder, __drem)
+strong_alias (__remainder_compat, __drem)
 compat_symbol (libm, __drem, dreml, GLIBC_2_0);
 #endif
--- a/sysdeps/m68k/m680x0/fpu/e_remainder.c
+++ b/sysdeps/m68k/m680x0/fpu/e_remainder.c
@@ -17,12 +17,26 @@
   <https://www.gnu.org/licenses/>.  */

 #include <math.h>
+#include <libm-alias-double.h>
 #include <libm-alias-finite.h>
 #include "mathimpl.h"
+#include "math_config.h"

 double
-__ieee754_remainder (double x, double y)
+__remainder (double x, double y)
 {
+  uint64_t hx = asuint64 (x);
+  uint64_t hy = asuint64 (y);
+
+  /* fmod(+-Inf,y) or fmod(x,0) */
+  if (__glibc_unlikely ((is_inf (hx) || y == 0.0f)
+			&& !is_nan (hy)
+			&& !is_nan (hx)))
+    return __math_invalid (x);
+
  return __m81_u(__ieee754_remainder)(x, y);
 }
+strong_alias (__remainder, __ieee754_remainder)
+versioned_symbol (libm, __remainder, remainder, GLIBC_2_43);
+libm_alias_double_other (__remainder, remainder)
 libm_alias_finite (__ieee754_remainder, __remainder)
--- a/sysdeps/mach/hurd/i386/libm.abilist
+++ b/sysdeps/mach/hurd/i386/libm.abilist
@@ -1328,6 +1328,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/aarch64/libm.abilist
+++ b/sysdeps/unix/sysv/linux/aarch64/libm.abilist
@@ -1294,6 +1294,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/alpha/libm.abilist
+++ b/sysdeps/unix/sysv/linux/alpha/libm.abilist
@@ -1453,6 +1453,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/arm/be/libm.abilist
+++ b/sysdeps/unix/sysv/linux/arm/be/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/arm/le/libm.abilist
+++ b/sysdeps/unix/sysv/linux/arm/le/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/hppa/libm.abilist
+++ b/sysdeps/unix/sysv/linux/hppa/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/i386/libm.abilist
+++ b/sysdeps/unix/sysv/linux/i386/libm.abilist
@@ -1335,6 +1335,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/m68k/coldfire/libm.abilist
+++ b/sysdeps/unix/sysv/linux/m68k/coldfire/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/m68k/m680x0/libm.abilist
+++ b/sysdeps/unix/sysv/linux/m68k/m680x0/libm.abilist
@@ -992,6 +992,7 @@ GLIBC_2.43 fmodf F
 GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/microblaze/be/libm.abilist
+++ b/sysdeps/unix/sysv/linux/microblaze/be/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/microblaze/le/libm.abilist
+++ b/sysdeps/unix/sysv/linux/microblaze/le/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/mips/mips32/libm.abilist
+++ b/sysdeps/unix/sysv/linux/mips/mips32/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/mips/mips64/libm.abilist
+++ b/sysdeps/unix/sysv/linux/mips/mips64/libm.abilist
@@ -1294,6 +1294,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/powerpc/powerpc32/fpu/libm.abilist
+++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/fpu/libm.abilist
@@ -1106,6 +1106,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/powerpc/powerpc32/nofpu/libm.abilist
+++ b/sysdeps/unix/sysv/linux/powerpc/powerpc32/nofpu/libm.abilist
@@ -1105,6 +1105,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/powerpc/powerpc64/be/libm.abilist
+++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/be/libm.abilist
@@ -1099,6 +1099,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/powerpc/powerpc64/le/libm.abilist
+++ b/sysdeps/unix/sysv/linux/powerpc/powerpc64/le/libm.abilist
@@ -1483,6 +1483,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/s390/s390-32/libm.abilist
+++ b/sysdeps/unix/sysv/linux/s390/s390-32/libm.abilist
@@ -1397,6 +1397,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/s390/s390-64/libm.abilist
+++ b/sysdeps/unix/sysv/linux/s390/s390-64/libm.abilist
@@ -1397,6 +1397,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/sh/be/libm.abilist
+++ b/sysdeps/unix/sysv/linux/sh/be/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/sh/le/libm.abilist
+++ b/sysdeps/unix/sysv/linux/sh/le/libm.abilist
@@ -959,6 +959,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/sparc/sparc32/libm.abilist
+++ b/sysdeps/unix/sysv/linux/sparc/sparc32/libm.abilist
@@ -1404,6 +1404,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/sparc/sparc64/libm.abilist
+++ b/sysdeps/unix/sysv/linux/sparc/sparc64/libm.abilist
@@ -1294,6 +1294,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/x86_64/64/libm.abilist
+++ b/sysdeps/unix/sysv/linux/x86_64/64/libm.abilist
@@ -1327,6 +1327,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F
--- a/sysdeps/unix/sysv/linux/x86_64/x32/libm.abilist
+++ b/sysdeps/unix/sysv/linux/x86_64/x32/libm.abilist
@@ -1327,6 +1327,7 @@ GLIBC_2.43 j0f F
 GLIBC_2.43 j1f F
 GLIBC_2.43 jnf F
 GLIBC_2.43 log10f F
+GLIBC_2.43 remainder F
 GLIBC_2.43 remainderf F
 GLIBC_2.43 y0f F
 GLIBC_2.43 y1f F