[PATCH v2] math: New generic fmaf implementation

[Adhemerval Zanella]

The current implementation relies on setting the rounding mode for
different calculations (FE_TOWARDZERO) to obtain correctly rounded
results. For most CPUs, this adds a significant performance overhead
since it requires executing a typically slow instruction (to get/set
the floating-point status), it necessitates flushing the pipeline,
and breaks some compiler assumptions/optimizations.

This patch introduces a new implementation originally written by
Szabolcs for musl, which removes the need to change rounding mode
and adds an extra step to round the numbers.

The original implementation adds some tests to handle underflow in
corner cases; however, it fails on the arm32 target. To fix it, I
kept the current implementation as a fallback to correctly raise
the expected exceptions.

I tested this implementation on various targets (x86_64, i686, arm,
aarch64, powerpc), including some by manually disabling the compiler
 instructions.

Performance-wise, it shows large improvements:

reciprocal-throughput       master       patched       improvement
x86_64 [1]                   58.09          7.96             7.33x
i686 [1]                    279.41         16.97            16.46x
aarch64 [2]                  26.09          4.10             6.35x
armhf [2]                    30.25          4.20             7.18x
powerpc [3]                   9.46          1.46             6.45x

latency                     master       patched       improvement
x86_64                       64.50         14.25             4.53x
i686                        304.39         61.04             4.99x
aarch64                      27.71          5.74             4.82x
armhf                        33.46          7.34             4.55x
powerpc                      10.96          2.65             4.13x

Checked on x86_64-linux-gnu and i686-linux-gnu with —disable-multi-arch,
and on arm-linux-gnueabihf.

[1] gcc 15.2.1, Zen3
[2] gcc 15.2.1, Neoverse N1
[3] gcc 15.2.1, POWER10

Signed-off-by: Szabolcs Nagy <nsz@gcc.gnu.org>
Co-authored-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
---
 sysdeps/i386/Makefile           |  1 +
 sysdeps/ieee754/dbl-64/s_fmaf.c | 58 +++++++++++++++------------------
 2 files changed, 28 insertions(+), 31 deletions(-)

diff --git a/sysdeps/i386/Makefile b/sysdeps/i386/Makefile
index 11ddbd402d..e53bb0c5cd 100644
--- a/sysdeps/i386/Makefile
+++ b/sysdeps/i386/Makefile
@@ -14,6 +14,7 @@ CFLAGS-s_erf.c += -fexcess-precision=standard
 CFLAGS-s_erfc.c += -fexcess-precision=standard
 CFLAGS-s_erf_common.c += -fexcess-precision=standard
 CFLAGS-s_fma.c += -fexcess-precision=standard
+CFLAGS-s_fmaf.c += -fexcess-precision=standard
 endif
 
 ifeq ($(subdir),gmon)
diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c
index 26ced15465..8add202b73 100644
--- a/sysdeps/ieee754/dbl-64/s_fmaf.c
+++ b/sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -19,16 +19,9 @@
 #define NO_MATH_REDIRECT
 #include <math.h>
 #include <fenv.h>
-#include <ieee754.h>
-#include <math-barriers.h>
-#include <fenv_private.h>
 #include <libm-alias-float.h>
-
-/* This implementation relies on double being more than twice as
-   precise as float and uses rounding to odd in order to avoid problems
-   with double rounding.
-   See a paper by Boldo and Melquiond:
-   http://www.lri.fr/~melquion/doc/08-tc.pdf  */
+#include <math-use-builtins.h>
+#include "math_config.h"
 
 float
 __fmaf (float x, float y, float z)
@@ -36,34 +29,37 @@ __fmaf (float x, float y, float z)
 #if USE_FMAF_BUILTIN
   return __builtin_fmaf (x, y, z);
 #else
-  /* Use generic implementation.  */
-  fenv_t env;
-
   /* Multiplication is always exact.  */
-  double temp = (double) x * (double) y;
+  double xy = (double) x * (double) y;
+  double result = xy + z;
 
-  /* Ensure correct sign of an exact zero result by performing the
-     addition in the original rounding mode in that case.  */
-  if (temp == (double) -z)
-    return (float) temp + z;
+  uint64_t u = asuint64 (result);
+  if (__glibc_likely ((u & 0xfffffff) != 0))
+    return result;
 
-  union ieee754_double u;
+  if ((u & 0x10000000) == 0
+      && ((u >> MANTISSA_WIDTH) & 0x7ff) > EXPONENT_BIAS - 126)
+    return result;
 
-  libc_feholdexcept_setround (&env, FE_TOWARDZERO);
+  if (result - xy == z && result - z == xy)
+    return result;
 
-  /* Perform addition with round to odd.  */
-  u.d = temp + (double) z;
-  /* Ensure the addition is not scheduled after fetestexcept call.  */
-  math_force_eval (u.d);
+  /*
+   * If result is inexact, and exactly halfway between two float values,
+   * we need to adjust the low-order bit in the direction of the error.
+   */
+  double err;
+  int neg = u >> 63;
+  if (neg == (z > xy))
+    err = xy - result + z;
+  else
+    err = z - result + xy;
 
-  /* Reset rounding mode and test for inexact simultaneously.  */
-  int j = libc_feupdateenv_test (&env, FE_INEXACT) != 0;
-
-  if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
-    u.ieee.mantissa1 |= j;
-
-  /* And finally truncation with round to nearest.  */
-  return (float) u.d;
+  if (neg == (err < 0))
+    u++;
+  else
+    u--;
+  return asdouble (u);
 #endif /* ! USE_FMAF_BUILTIN  */
 }
 #ifndef __fmaf
-- 
2.43.0