Core: Resolve EIGEN_SCALAR_MADD_USE_FMA after EIGEN_VECTORIZE_FMA is known

libeigen/eigen!2992

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/Core/util/ConfigureVectorization.h

diff --git a/Eigen/src/Core/util/ConfigureVectorization.h b/Eigen/src/Core/util/ConfigureVectorization.h
index 3f7612e..07e682c 100644
--- a/Eigen/src/Core/util/ConfigureVectorization.h
+++ b/Eigen/src/Core/util/ConfigureVectorization.h
@@ -664,6 +664,20 @@
 // IWYU pragma: private
 #include "../InternalHeaderCheck.h"
 
+/** Whether numext::madd uses std::fma for scalars. Defaults to the hardware: fused where a single
+ * instruction exists, which keeps the scalar and vectorized paths consistent, and unfused otherwise,
+ * where software fma costs 2-3x on Intel and up to 30x on WASM. Resolved here rather than in Macros.h,
+ * which Eigen/Core includes first, because the architecture branches above settle
+ * EIGEN_VECTORIZE_FMA -- the ARM one only a few lines up.
+ */
+#ifndef EIGEN_SCALAR_MADD_USE_FMA
+#ifdef EIGEN_VECTORIZE_FMA
+#define EIGEN_SCALAR_MADD_USE_FMA 1
+#else
+#define EIGEN_SCALAR_MADD_USE_FMA 0
+#endif
+#endif
+
 namespace Eigen {
 
 inline static const char* SimdInstructionSetsInUse(void) {

Eigen/src/Core/util/Macros.h

diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 8e6e0e7..35fdbaf 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -61,26 +61,6 @@
 #define EIGEN_STACK_ALLOCATION_LIMIT_WAS_DEFAULTED
 #endif
 
-/* Specify whether to use std::fma for scalar multiply-add instructions.
- *
- * On machines that have FMA as a single instruction, this will generally
- * improve precision without significant performance implications.
- *
- * Without a single instruction, performance has been found to be reduced 2-3x
- * on Intel CPUs, and up to 30x for WASM.
- *
- * If unspecified, defaults to using FMA if hardware support is available.
- * The default should be used in most cases to ensure consistency between
- * vectorized and non-vectorized paths.
- */
-#ifndef EIGEN_SCALAR_MADD_USE_FMA
-#ifdef EIGEN_VECTORIZE_FMA
-#define EIGEN_SCALAR_MADD_USE_FMA 1
-#else
-#define EIGEN_SCALAR_MADD_USE_FMA 0
-#endif
-#endif
-
 //------------------------------------------------------------------------------------------
 // Compiler identification, EIGEN_COMP_*
 //------------------------------------------------------------------------------------------

test/numext.cpp

diff --git a/test/numext.cpp b/test/numext.cpp
index 94b4f03..50db67b 100644
--- a/test/numext.cpp
+++ b/test/numext.cpp
@@ -261,6 +261,24 @@
   }
 }
 
+// internal::twoprod splits x*y into hi + lo with x*y = hi + lo exactly, which holds only if its
+// multiply-add is genuinely fused: a non-fused x*y - hi folds against the already-rounded product and
+// silently yields lo = 0, leaving every double-word computation built on it with no low word at all.
+template <typename T>
+void check_twoprod() {
+  for (int k = 0; k < 100; ++k) {
+    const T x = internal::random<T>(T(-1), T(1));
+    const T y = internal::random<T>(T(-1), T(1));
+    T hi, lo;
+    internal::twoprod(x, y, hi, lo);
+    VERIFY_IS_EQUAL(hi, x * y);
+    // fma is exact by IEEE-754 contract, so this is the definition of the low word rather than an
+    // independent approximation of it.
+    EIGEN_USING_STD(fma);
+    VERIFY_IS_EQUAL(lo, fma(x, y, -hi));
+  }
+}
+
 template <typename T>
 void check_arg() {
   typedef typename NumTraits<T>::Real Real;
@@ -600,6 +618,9 @@
     CALL_SUBTEST(check_complex_sign<std::complex<float>>());
     CALL_SUBTEST(check_complex_sign<std::complex<double>>());
 
+    CALL_SUBTEST(check_twoprod<float>());
+    CALL_SUBTEST(check_twoprod<double>());
+
     CALL_SUBTEST(check_arg<std::complex<float>>());
     CALL_SUBTEST(check_arg<std::complex<double>>());