math: Remove the error handling wrapper from fmod and fmodf

The error handling is moved to sysdeps/ieee754 version with no SVID
support.  The compatibility symbol versions still use the wrapper
with SVID error handling around the new code.  There is no new symbol
version nor compatibility code on !LIBM_SVID_COMPAT targets
(e.g. riscv).

The ia64 is unchanged, since it still uses the arch specific
__libm_error_region on its implementation.  For both i686 and m68k,
which provive arch specific implementation, wrappers are added so
no new symbol are added (which would require to change the
implementations).

It shows an small improvement, the results for fmod:

  Architecture     | Input           | master   | patch
  -----------------|-----------------|----------|--------
  x86_64 (Ryzen 9) | subnormals      | 12.5049  | 9.40992
  x86_64 (Ryzen 9) | normal          | 296.939  | 296.738
  x86_64 (Ryzen 9) | close-exponents | 16.0244  | 13.119
  aarch64 (N1)     | subnormal       | 6.81778  | 4.33313
  aarch64 (N1)     | normal          | 155.620  | 152.915
  aarch64 (N1)     | close-exponents | 8.21306  | 5.76138
  armhf (N1)       | subnormal       | 15.1083  | 14.5746
  armhf (N1)       | normal          | 244.833  | 241.738
  armhf (N1)       | close-exponents | 21.8182  | 22.457

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

8 files changed, 53 insertions, 7 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c
index e661ca1ff8..caae4e47e2 100644
--- a/sysdeps/ieee754/dbl-64/e_fmod.c
+++ b/sysdeps/ieee754/dbl-64/e_fmod.c
@@ -16,7 +16,9 @@
    License along with the GNU C Library; if not, see
    <https://www.gnu.org/licenses/>.  */
 
+#include <libm-alias-double.h>
 #include <libm-alias-finite.h>
+#include <math-svid-compat.h>
 #include <math.h>
 #include "math_config.h"
 
@@ -55,7 +57,7 @@
      }  */
 
 double
-__ieee754_fmod (double x, double y)
+__fmod (double x, double y)
 {
   uint64_t hx = asuint64 (x);
   uint64_t hy = asuint64 (y);
@@ -67,11 +69,16 @@ __ieee754_fmod (double x, double y)
 
   /* Special cases:
      - If x or y is a Nan, NaN is returned.
-     - If x is an inifinity, a NaN is returned.
+     - If x is an inifinity, a NaN is returned and EDOM is set.
      - If y is zero, Nan is returned.
      - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
-  if (__glibc_unlikely (hy == 0	|| hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
-    return (x * y) / (x * y);
+  if (__glibc_unlikely (hy == 0
+			|| hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
+    {
+      if (is_nan (hx) || is_nan (hy))
+	return (x * y) / (x * y);
+      return __math_edom ((x * y) / (x * y));
+    }
 
   if (__glibc_unlikely (hx <= hy))
     {
@@ -153,4 +160,11 @@ __ieee754_fmod (double x, double y)
 
   return make_double (mx, ey, sx);
 }
+strong_alias (__fmod, __ieee754_fmod)
 libm_alias_finite (__ieee754_fmod, __fmod)
+#if LIBM_SVID_COMPAT
+versioned_symbol (libm, __fmod, fmod, GLIBC_2_38);
+libm_alias_double_other (__fmod, fmod)
+#else
+libm_alias_double (__fmod, fmod)
+#endif
diff --git a/sysdeps/ieee754/dbl-64/math_config.h b/sysdeps/ieee754/dbl-64/math_config.h
index 2049cea3f7..499ca7591b 100644
--- a/sysdeps/ieee754/dbl-64/math_config.h
+++ b/sysdeps/ieee754/dbl-64/math_config.h
@@ -170,6 +170,9 @@ attribute_hidden double __math_invalid (double);
 
 /* Error handling using output checking, only for errno setting.  */
 
+/* Check if the result generated a demain error.  */
+attribute_hidden double __math_edom (double x);
+
 /* Check if the result overflowed to infinity.  */
 attribute_hidden double __math_check_oflow (double);
 /* Check if the result underflowed to 0.  */
diff --git a/sysdeps/ieee754/dbl-64/math_err.c b/sysdeps/ieee754/dbl-64/math_err.c
index 5f5f965d3e..377a4caad2 100644
--- a/sysdeps/ieee754/dbl-64/math_err.c
+++ b/sysdeps/ieee754/dbl-64/math_err.c
@@ -33,6 +33,12 @@ with_errno (double y, int e)
 #define with_errno(x, e) (x)
 #endif
 
+attribute_hidden double
+__math_edom (double y)
+{
+  return with_errno (y, EDOM);
+}
+
 /* NOINLINE reduces code size.  */
 NOINLINE static double
 xflow (uint32_t sign, double y)
diff --git a/sysdeps/ieee754/dbl-64/w_fmod.c b/sysdeps/ieee754/dbl-64/w_fmod.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/w_fmod.c
@@ -0,0 +1 @@
+/* Not needed.  */
diff --git a/sysdeps/ieee754/flt-32/e_fmodf.c b/sysdeps/ieee754/flt-32/e_fmodf.c
index 4482b8c0ac..763900efda 100644
--- a/sysdeps/ieee754/flt-32/e_fmodf.c
+++ b/sysdeps/ieee754/flt-32/e_fmodf.c
@@ -17,6 +17,8 @@
    <https://www.gnu.org/licenses/>.  */
 
 #include <libm-alias-finite.h>
+#include <libm-alias-float.h>
+#include <math-svid-compat.h>
 #include <math.h>
 #include "math_config.h"
 
@@ -55,7 +57,7 @@
      }  */
 
 float
-__ieee754_fmodf (float x, float y)
+__fmodf (float x, float y)
 {
   uint32_t hx = asuint (x);
   uint32_t hy = asuint (y);
@@ -70,8 +72,13 @@ __ieee754_fmodf (float x, float y)
      - If x is an inifinity, a NaN is returned.
      - If y is zero, Nan is returned.
      - If x is +0/-0, and y is not zero, +0/-0 is returned.  */
-  if (__glibc_unlikely (hy == 0	|| hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
-    return (x * y) / (x * y);
+  if (__glibc_unlikely (hy == 0
+			|| hx >= EXPONENT_MASK || hy > EXPONENT_MASK))
+    {
+      if (is_nan (hx) || is_nan (hy))
+	return (x * y) / (x * y);
+      return __math_edomf ((x * y) / (x * y));
+    }
 
   if (__glibc_unlikely (hx <= hy))
     {
@@ -152,4 +159,11 @@ __ieee754_fmodf (float x, float y)
 
   return make_float (mx, ey, sx);
 }
+strong_alias (__fmodf, __ieee754_fmodf)
+#if LIBM_SVID_COMPAT
+versioned_symbol (libm, __fmodf, fmodf, GLIBC_2_38);
+libm_alias_float_other (__fmod, fmod)
+#else
+libm_alias_float (__fmod, fmod)
+#endif
 libm_alias_finite (__ieee754_fmodf, __fmodf)
diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
index 829430ea28..a7fb332767 100644
--- a/sysdeps/ieee754/flt-32/math_config.h
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -158,6 +158,7 @@ attribute_hidden float __math_uflowf (uint32_t);
 attribute_hidden float __math_may_uflowf (uint32_t);
 attribute_hidden float __math_divzerof (uint32_t);
 attribute_hidden float __math_invalidf (float);
+attribute_hidden float __math_edomf (float x);
 
 /* Shared between expf, exp2f, exp10f, and powf.  */
 #define EXP2F_TABLE_BITS 5
diff --git a/sysdeps/ieee754/flt-32/math_errf.c b/sysdeps/ieee754/flt-32/math_errf.c
index 0389ce0a81..7245ab8d96 100644
--- a/sysdeps/ieee754/flt-32/math_errf.c
+++ b/sysdeps/ieee754/flt-32/math_errf.c
@@ -31,6 +31,12 @@ with_errnof (float y, int e)
 # define with_errnof(x, e) (x)
 #endif
 
+attribute_hidden float
+__math_edomf (float y)
+{
+  return with_errnof (y, EDOM);
+}
+
 /* NOINLINE prevents fenv semantics breaking optimizations.  */
 NOINLINE static float
 xflowf (uint32_t sign, float y)
diff --git a/sysdeps/ieee754/flt-32/w_fmodf.c b/sysdeps/ieee754/flt-32/w_fmodf.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/w_fmodf.c
@@ -0,0 +1 @@
+/* Not needed.  */