[PATCH 5/5] math: New generic fmaf implementation
Adhemerval Zanella
adhemerval.zanella@linaro.org
Thu Nov 13 12:58:22 GMT 2025
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 8.78 6.62x
i686 [1] 279.41 18.30 15.26x
aarch64 [2] 26.09 4.24 6.16x
armhf [2] 30.25 4.93 6.13x
powerpc [3] 9.46 1.87 5.04x
latency master patched improvement
x86_64 64.50 14.35 4.49x
i686 304.39 64.43 4.72x
aarch64 27.71 5.91 4.69x
armhf 33.46 7.26 4.60x
powerpc 10.96 3.43 3.19x
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 | 87 +++++++++++++++++++++++----------
2 files changed, 61 insertions(+), 27 deletions(-)
diff --git a/sysdeps/i386/Makefile b/sysdeps/i386/Makefile
index 11ddbd402d0..e53bb0c5cd4 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 7bf9941dd88..8861d18c214 100644
--- a/sysdeps/ieee754/dbl-64/s_fmaf.c
+++ b/sysdeps/ieee754/dbl-64/s_fmaf.c
@@ -19,11 +19,14 @@
#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>
#include <math-use-builtins.h>
+#include "math_config.h"
+
+#if !USE_FMAF_BUILTIN
+# include <ieee754.h>
+# include <math-barriers.h>
+# include <fenv_private.h>
/* This implementation relies on double being more than twice as
precise as float and uses rounding to odd in order to avoid problems
@@ -31,40 +34,70 @@
See a paper by Boldo and Melquiond:
http://www.lri.fr/~melquion/doc/08-tc.pdf */
+static __attribute_noinline__ float
+fmaf_fallback (double xy, double z)
+{
+ fenv_t env;
+ union ieee754_double u;
+
+ libc_feholdexcept_setround (&env, FE_TOWARDZERO);
+
+ /* Perform addition with round to odd. */
+ u.d = xy + (double) z;
+ /* Ensure the addition is not scheduled after fetestexcept call. */
+ math_force_eval (u.d);
+
+ /* 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 u.d;
+}
+#endif
+
float
__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);
+ /* Common case: The double precision result is fine. */
+ if ((u & 0x1fffffff) != 0x10000000 || /* not a halfway case */
+ (result - xy == z && result - z == xy) || /* exact */
+ __fegetround () != FE_TONEAREST) /* not round-to-nearest */
+ {
+ /* Underflow may not be raised correctly, example:
+ fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f) */
+ int e = u >> MANTISSA_WIDTH & 0x7ff;
+ if (__glibc_unlikely (e <= EXPONENT_BIAS - 126
+ && e >= EXPONENT_BIAS - 149))
+ return fmaf_fallback (xy, z);
+ return result;
+ }
- union ieee754_double u;
+ /*
+ * 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;
- libc_feholdexcept_setround (&env, FE_TOWARDZERO);
-
- /* 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);
-
- /* 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
More information about the Libc-alpha
mailing list