Core: Preserve double-word products with scalar FMA disabled libeigen/eigen!3005 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h b/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h index 27c0239..7840951 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h
@@ -35,6 +35,18 @@ } #ifdef EIGEN_VECTORIZE_FMA +// Given x, y, and xy = fl(x*y), return the residual such that x*y = xy + residual exactly. +template <typename Packet, std::enable_if_t<!is_scalar<Packet>::value, int> = 0> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) { + return pmsub(x, y, xy); +} + +template <typename Scalar, std::enable_if_t<is_scalar<Scalar>::value, int> = 0> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low(const Scalar& x, const Scalar& y, const Scalar& xy) { + // Error-free products require FMA even when EIGEN_SCALAR_MADD_USE_FMA disables fusion in scalar madd. + return numext::fma(x, y, Scalar(-xy)); +} + // This function implements the extended precision product of // a pair of floating point numbers. Given {x, y}, it computes the pair // {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and @@ -42,14 +54,7 @@ template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) { p_hi = pmul(x, y); - p_lo = pmsub(x, y, p_hi); -} - -// A version of twoprod that takes x, y, and fl(x*y) as input and returns the p_lo such that -// x * y = xy + p_lo holds exactly. -template <typename Packet> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) { - return pmsub(x, y, xy); + p_lo = twoprod_low(x, y, p_hi); } #else
diff --git a/Eigen/src/Core/util/ConfigureVectorization.h b/Eigen/src/Core/util/ConfigureVectorization.h index 07e682c..084960b 100644 --- a/Eigen/src/Core/util/ConfigureVectorization.h +++ b/Eigen/src/Core/util/ConfigureVectorization.h
@@ -666,9 +666,10 @@ /** 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. + * where software fma costs 2-3x on Intel and up to 30x on WASM. Error-free transformations such as + * internal::twoprod are exact only when fused, so they call numext::fma regardless of this setting. + * 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
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index da6f25b..6a4b588 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -274,6 +274,7 @@ ei_add_test(lru_cache) ei_add_test(maxsizevector) ei_add_test(numext) +ei_add_test(numext_fma_off) ei_add_test(sizeof) ei_add_test(dynalloc) ei_add_test(nomalloc)
diff --git a/test/numext.cpp b/test/numext.cpp index 50db67b..5375652 100644 --- a/test/numext.cpp +++ b/test/numext.cpp
@@ -261,24 +261,66 @@ } } -// 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. +// The FMA implementation must retain the product's rounding error even when scalar madd is unfused. +template <typename T> +void check_twoprod_pair(const T& x, const T& y) { + 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)); + VERIFY_IS_EQUAL(internal::twoprod_low(x, y, hi), fma(x, y, -hi)); +} + template <typename T> void check_twoprod() { + // (1 + epsilon) * (1 - epsilon) = 1 - epsilon^2 exposes a lost low word in every IEEE type. + const T epsilon = NumTraits<T>::epsilon(); + check_twoprod_pair(T(1) + epsilon, T(1) - epsilon); 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)); + check_twoprod_pair(internal::random<T>(T(-1), T(1)), internal::random<T>(T(-1), T(1))); } } +#ifdef EIGEN_VECTORIZE_FMA +// Custom scalars can return a lazy expression from unary minus. +struct TwoprodScalar { + struct Negation { + const TwoprodScalar& operand; + }; + + TwoprodScalar() = default; + explicit TwoprodScalar(double v) : value(v) {} + TwoprodScalar(const Negation& negation) : value(-negation.operand.value) {} + + Negation operator-() const { return {*this}; } + TwoprodScalar operator*(const TwoprodScalar& other) const { return TwoprodScalar(value * other.value); } + + friend TwoprodScalar fma(const TwoprodScalar& x, const TwoprodScalar& y, const TwoprodScalar& z) { + return TwoprodScalar(std::fma(x.value, y.value, z.value)); + } + + double value = 0; +}; + +void check_twoprod_negation_expression() { + STATIC_CHECK(internal::is_scalar<TwoprodScalar>::value); + STATIC_CHECK(internal::has_fma<TwoprodScalar>::value); + STATIC_CHECK((!std::is_same<decltype(-std::declval<TwoprodScalar>()), TwoprodScalar>::value)); + const double epsilon = NumTraits<double>::epsilon(); + for (double sign : {-1.0, 1.0}) { + const TwoprodScalar x(sign * (1 + epsilon)), y(1 - epsilon); + TwoprodScalar hi, lo; + internal::twoprod(x, y, hi, lo); + VERIFY_IS_EQUAL(hi.value, sign); + VERIFY_IS_EQUAL(lo.value, -sign * epsilon * epsilon); + VERIFY_IS_EQUAL(internal::twoprod_low(x, y, hi).value, -sign * epsilon * epsilon); + } +} +#endif + template <typename T> void check_arg() { typedef typename NumTraits<T>::Real Real; @@ -620,6 +662,12 @@ CALL_SUBTEST(check_twoprod<float>()); CALL_SUBTEST(check_twoprod<double>()); + CALL_SUBTEST(check_twoprod<long double>()); + CALL_SUBTEST(check_twoprod<half>()); + CALL_SUBTEST(check_twoprod<bfloat16>()); +#ifdef EIGEN_VECTORIZE_FMA + CALL_SUBTEST(check_twoprod_negation_expression()); +#endif CALL_SUBTEST(check_arg<std::complex<float>>()); CALL_SUBTEST(check_arg<std::complex<double>>());
diff --git a/test/numext_fma_off.cpp b/test/numext_fma_off.cpp new file mode 100644 index 0000000..0619b94 --- /dev/null +++ b/test/numext_fma_off.cpp
@@ -0,0 +1,6 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#undef EIGEN_SCALAR_MADD_USE_FMA +#define EIGEN_SCALAR_MADD_USE_FMA 0 +#include "numext.cpp" // NOLINT(bugprone-suspicious-include): Compile the suite with a different configuration.