From c4bed222b2afe6fef76b0a0fb8e88acff9cc32d5 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com> Date: Wed, 8 Jul 2026 22:59:20 -0700 Subject: [PATCH] Core: build packet sign masks from integer bits to survive fast-math libeigen/eigen!2699 Co-authored-by: Rasmus Munk Larsen --- Eigen/src/Core/GenericPacketMath.h | 25 +++++++++- .../arch/Default/GenericPacketMathComplex.h | 5 +- .../arch/Default/GenericPacketMathFunctions.h | 2 +- .../Core/arch/Default/GenericPacketMathPow.h | 4 +- .../Core/arch/Default/GenericPacketMathTrig.h | 20 ++++---- Eigen/src/Core/arch/SSE/Complex.h | 4 +- Eigen/src/LU/arch/InverseSize4.h | 10 ++-- test/fastmath.cpp | 48 +++++++++++++++++++ .../SpecialFunctions/SpecialFunctionsImpl.h | 3 +- 9 files changed, 96 insertions(+), 25 deletions(-) diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index 89e3adea2..cba85dce8 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h @@ -819,7 +819,28 @@ EIGEN_DEVICE_FUNC inline Packet pset1(const typename unpacket_traits::ty /** \internal \returns a packet with constant coefficients set from bits */ template -EIGEN_DEVICE_FUNC inline Packet pset1frombits(BitsType a); +EIGEN_DEVICE_FUNC inline Packet pset1frombits(BitsType a) { + using Scalar = typename unpacket_traits::type; + return pset1(numext::bit_cast(a)); +} + +/** \internal \returns a packet with all coefficients set to -0.0, i.e. with only the sign bit set. + * + * The mask is deliberately constructed from the integer sign-bit pattern via pset1frombits + * instead of the floating-point literal -Scalar(0): under fast-math flags (-ffast-math implies + * -fno-signed-zeros) compilers may treat -0.0 and +0.0 as interchangeable, and e.g. GCC's + * value numbering substitutes a splat of -0.0 with a nearby splat of +0.0, silently zeroing + * the mask and corrupting sign manipulation of non-zero values. Architectures that specialize + * pset1frombits keep the constant in the integer domain, where no floating-point + * simplification applies. See https://gitlab.com/libeigen/eigen/-/merge_requests/2698. + */ +template +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet psignmask() { + using Scalar = typename unpacket_traits::type; + using Bits = typename numext::get_integer_by_size::unsigned_type; + constexpr Bits kSignBit = static_cast(Bits(1) << (CHAR_BIT * sizeof(Scalar) - 1)); + return pset1frombits(kSignBit); +} template ::value, int> = 0> EIGEN_DEVICE_FUNC inline Scalar pload1_scalar(const Scalar* a) { @@ -1589,7 +1610,7 @@ EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet patan2(const Packet& y, const Packe // See https://en.cppreference.com/w/cpp/numeric/math/atan2 // for how corner cases are supposed to be handled according to the // IEEE floating-point standard (IEC 60559). - const Packet kSignMask = pset1(-Scalar(0)); + const Packet kSignMask = psignmask(); const Packet kZero = pzero(x); const Packet kOne = pset1(Scalar(1)); const Packet kPi = pset1(Scalar(EIGEN_PI)); diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h index 7c0a3ea7c..94657e9ef 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h @@ -121,7 +121,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pexp_complex(const Pa // prevent inf * 0 = NaN. The vectorized sincos may compute exact zero // for near-zero values like cos(pi/2), and inf * +-1 = +-inf is correct. // The y=0 case is handled separately below. - RealPacket cisy_sign_one = por(pand(cisy, pset1(RealScalar(-0.0))), pset1(RealScalar(1))); + RealPacket cisy_sign_one = por(pand(cisy, psignmask()), pset1(RealScalar(1))); RealPacket expx_inf_y_finite = pand(pcmp_eq(expx, cst_pos_inf), pcmp_lt(pabs(y), cst_pos_inf)); cisy = pselect(expx_inf_y_finite, cisy_sign_one, cisy); @@ -205,7 +205,8 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psqrt_complex(const P // Step 4. Compute solution for inputs with negative real part: // [|eta0|, sign(y0)*rho0, |eta1|, sign(y1)*rho1] - const RealPacket cst_imag_sign_mask = pset1(Scalar(RealScalar(0.0), RealScalar(-0.0))).v; + // [+0.0, -0.0, ...]: the sign bit of the imaginary (odd) lanes only. + const RealPacket cst_imag_sign_mask = pandnot(psignmask(), real_mask); RealPacket imag_signs = pand(a.v, cst_imag_sign_mask); Packet negative_real_result; // Notice that rho is positive, so taking its absolute value is a noop. diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index 4fc5c8c60..0463717c7 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h @@ -777,7 +777,7 @@ template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet generic_ceil(const Packet& a) { using Scalar = typename unpacket_traits::type; const Packet cst_1 = pset1(Scalar(1)); - const Packet sign_mask = pset1(static_cast(-0.0)); + const Packet sign_mask = psignmask(); Packet rint_a = generic_rint(a); // if rint(a) < a, then rint(a) == floor(a) Packet mask = pcmp_lt(rint_a, a); diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathPow.h b/Eigen/src/Core/arch/Default/GenericPacketMathPow.h index 2642e29b4..2b3188a90 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathPow.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathPow.h @@ -58,10 +58,8 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet cbrt_decompose(const template EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet cbrt_special_cases_and_sign(const Packet& x, const Packet& abs_root) { - typedef typename unpacket_traits::type Scalar; - // Set sign. - const Packet sign_mask = pset1(Scalar(-0.0)); + const Packet sign_mask = psignmask(); const Packet x_sign = pand(sign_mask, x); Packet root = por(x_sign, abs_root); diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h index dd6826445..c15a013ab 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h @@ -96,7 +96,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS const Packet cst_2oPI = pset1(0.636619746685028076171875f); // 2/PI const Packet cst_rounding_magic = pset1(12582912); // 2^23 for rounding const PacketI csti_1 = pset1(1); - const Packet cst_sign_mask = pset1frombits(static_cast(0x80000000u)); + const Packet cst_sign_mask = psignmask(); Packet x = pabs(_x); @@ -294,7 +294,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS typedef typename unpacket_traits::integer_packet PacketI; typedef typename unpacket_traits::type ScalarI; - const Packet cst_sign_mask = pset1frombits(static_cast(0x8000000000000000u)); + const Packet cst_sign_mask = psignmask(); // If the argument is smaller than this value, use a simpler argument reduction const double small_th = 15; @@ -579,7 +579,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet generic_atan(const Pa constexpr Scalar kPiOverTwo = static_cast(EIGEN_PI / 2); - const Packet cst_signmask = pset1(Scalar(-0.0)); + const Packet cst_signmask = psignmask(); const Packet cst_one = pset1(Scalar(1)); const Packet cst_pi_over_two = pset1(kPiOverTwo); @@ -679,7 +679,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS T ptanh_float(const T& x) { const T q = ppolevl::run(x2, beta); const T small_tanh = pmadd(x3, pdiv(p, q), x); - const T sign_mask = pset1(-0.0f); + const T sign_mask = psignmask(); const T abs_x = pandnot(x, sign_mask); constexpr float kSmallThreshold = 1.25f; const T large_mask = pcmp_lt(pset1(kSmallThreshold), abs_x); @@ -790,7 +790,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_float(const Pa const Packet x_eq_one = pcmp_eq(one, pabs(x)); const Packet x_gt_one = pcmp_lt(one, pabs(x)); - const Packet sign_mask = pset1(-0.0f); + const Packet sign_mask = psignmask(); const Packet x_sign = pand(sign_mask, x); const Packet inf = pset1(std::numeric_limits::infinity()); return por(x_gt_one, pselect(x_eq_one, por(x_sign, inf), pselect(x_gt_half, r, p))); @@ -828,7 +828,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet patanh_double(const P const Packet x_eq_one = pcmp_eq(one, pabs(x)); const Packet x_gt_one = pcmp_lt(one, pabs(x)); - const Packet sign_mask = pset1(-0.0); + const Packet sign_mask = psignmask(); const Packet x_sign = pand(sign_mask, x); const Packet inf = pset1(std::numeric_limits::infinity()); return por(x_gt_one, pselect(x_eq_one, por(x_sign, inf), pselect(x_gt_half, y_large, y_small))); @@ -848,7 +848,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_float(const Pac typedef typename unpacket_traits::type Scalar; static_assert(std::is_same::value, "Scalar type must be float"); - const Packet sign_mask = pset1(-0.0f); + const Packet sign_mask = psignmask(); const Packet abs_x = pandnot(x, sign_mask); const Packet x_sign = pand(x, sign_mask); @@ -888,7 +888,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet psinh_double(const Pa typedef typename unpacket_traits::type Scalar; static_assert(std::is_same::value, "Scalar type must be double"); - const Packet sign_mask = pset1(-0.0); + const Packet sign_mask = psignmask(); const Packet abs_x = pandnot(x, sign_mask); const Packet x_sign = pand(x, sign_mask); @@ -990,7 +990,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pcosh_double(const Pa */ template EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Packet& x) { - const Packet sign_mask = pset1(-0.0f); + const Packet sign_mask = psignmask(); const Packet abs_x = pandnot(x, sign_mask); const Packet x_sign = pand(x, sign_mask); const Packet one = pset1(1.0f); @@ -1015,7 +1015,7 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_float(const Pa template EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pasinh_double(const Packet& x) { - const Packet sign_mask = pset1(-0.0); + const Packet sign_mask = psignmask(); const Packet abs_x = pandnot(x, sign_mask); const Packet x_sign = pand(x, sign_mask); const Packet one = pset1(1.0); diff --git a/Eigen/src/Core/arch/SSE/Complex.h b/Eigen/src/Core/arch/SSE/Complex.h index da88b3b2a..47bcf532a 100644 --- a/Eigen/src/Core/arch/SSE/Complex.h +++ b/Eigen/src/Core/arch/SSE/Complex.h @@ -105,7 +105,7 @@ EIGEN_STRONG_INLINE Packet2cf pmul(const Packet2cf& a, const Packet2cf& b) { #ifdef EIGEN_VECTORIZE_SSE3 __m128 result = _mm_addsub_ps(_mm_mul_ps(tmp2, b.v), tmp1); #else - const __m128 mask = _mm_setr_ps(-0.0f, 0.0f, -0.0f, 0.0f); + const __m128 mask = _mm_castsi128_ps(_mm_setr_epi32(0x80000000, 0x00000000, 0x80000000, 0x00000000)); __m128 result = _mm_add_ps(_mm_mul_ps(tmp2, b.v), _mm_xor_ps(tmp1, mask)); #endif #endif @@ -299,7 +299,7 @@ EIGEN_STRONG_INLINE Packet1cd pmul(const Packet1cd& a, const Packet1cd& b) { #ifdef EIGEN_VECTORIZE_SSE3 __m128d result = _mm_addsub_pd(_mm_mul_pd(tmp2, b.v), tmp1); #else - const __m128d mask = _mm_setr_pd(-0.0, 0.0); + const __m128d mask = _mm_castsi128_pd(_mm_set_epi64x(0x0, 0x8000000000000000ull)); __m128d result = _mm_add_pd(_mm_mul_pd(tmp2, b.v), _mm_xor_pd(tmp1, mask)); #endif #endif diff --git a/Eigen/src/LU/arch/InverseSize4.h b/Eigen/src/LU/arch/InverseSize4.h index 38bf0cc8a..76e78db4a 100644 --- a/Eigen/src/LU/arch/InverseSize4.h +++ b/Eigen/src/LU/arch/InverseSize4.h @@ -148,7 +148,9 @@ struct compute_inverse_size4(0x80000000u); + EIGEN_ALIGN_MAX const float sign_mask[4] = {0.0f, neg_zero, neg_zero, 0.0f}; const Packet4f p4f_sign_PNNP = pload(sign_mask); rd = pxor(rd, p4f_sign_PNNP); iA = pmul(iA, rd); @@ -324,8 +326,10 @@ struct compute_inverse_size4(0x8000000000000000ull); + EIGEN_ALIGN_MAX const double sign_mask1[2] = {0.0, neg_zero}; + EIGEN_ALIGN_MAX const double sign_mask2[2] = {neg_zero, 0.0}; const Packet2d sign_PN = pload(sign_mask1); const Packet2d sign_NP = pload(sign_mask2); d1 = pxor(rd, sign_PN); diff --git a/test/fastmath.cpp b/test/fastmath.cpp index 6e9fcfcd7..0db99b36f 100644 --- a/test/fastmath.cpp +++ b/test/fastmath.cpp @@ -287,6 +287,52 @@ void check_complex_fastmath() { check_complex_householder_qr(); } +// The packet implementations of these functions manipulate signs of non-zero values through a +// -0.0 bitmask. Under fast-math flags compilers consider -0.0 and +0.0 interchangeable and may +// substitute one such constant for the other (GCC's value numbering does this on RISC-V), which +// silently zeroes the mask unless it is constructed from integer bits (see psignmask and +// https://gitlab.com/libeigen/eigen/-/merge_requests/2698). Sign handling of *non-zero* inputs +// and outputs is not relaxed by fast-math, so these checks must hold. +template +void check_sign_dependent_functions() { + typedef Array ArrayType; + const Index n = 64; + + ArrayType ya(n), xa(n); + for (Index i = 0; i < n; ++i) { + Scalar s = Scalar(1) + Scalar(i) / Scalar(n); + // All four quadrants. + ya[i] = (i & 1) ? s : -s; + xa[i] = (i & 2) ? Scalar(2) * s : Scalar(-2) * s; + } + // Cover the |x| == |y| special path of patan2 in all four quadrants. + ya[0] = Scalar(1), xa[0] = Scalar(1); + ya[1] = Scalar(1), xa[1] = Scalar(-1); + ya[2] = Scalar(-1), xa[2] = Scalar(1); + ya[3] = Scalar(-1), xa[3] = Scalar(-1); + + const ArrayType atan2_result = ya.atan2(xa); + for (Index i = 0; i < n; ++i) { + VERIFY_IS_APPROX(atan2_result[i], std::atan2(ya[i], xa[i])); + } + + // Mixed-sign inputs, no exact zeros, |w| up to ~3 to hit the small- and large-|x| branches. + // atanh and cbrt also use the sign mask but are not checked here: the atanh mask only + // affects the |x| == 1 -> +/-inf path, which is meaningless under -ffinite-math-only, and + // vectorized cbrt is currently broken under clang fast-math for unrelated reasons. + const ArrayType w = ArrayType::LinSpaced(n, Scalar(-3), Scalar(3)) + Scalar(0.017); + const ArrayType atan_result = w.atan(); + const ArrayType sinh_result = w.sinh(); + const ArrayType tanh_result = w.tanh(); + const ArrayType asinh_result = w.asinh(); + for (Index i = 0; i < n; ++i) { + VERIFY_IS_APPROX(atan_result[i], std::atan(w[i])); + VERIFY_IS_APPROX(sinh_result[i], std::sinh(w[i])); + VERIFY_IS_APPROX(tanh_result[i], std::tanh(w[i])); + VERIFY_IS_APPROX(asinh_result[i], std::asinh(w[i])); + } +} + EIGEN_DECLARE_TEST(fastmath) { std::cout << "*** float *** \n\n"; check_inf_nan(true); @@ -301,4 +347,6 @@ EIGEN_DECLARE_TEST(fastmath) { CALL_SUBTEST_1(check_complex_fastmath()); CALL_SUBTEST_2(check_complex_fastmath()); + CALL_SUBTEST_3(check_sign_dependent_functions()); + CALL_SUBTEST_4(check_sign_dependent_functions()); } diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index 10351ce7a..e7e9296d1 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -612,8 +612,7 @@ struct erf_impl { template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T flipsign(const T& should_flipsign, const T& x) { - typedef typename unpacket_traits::type Scalar; - const T sign_mask = pset1(Scalar(-0.0)); + const T sign_mask = psignmask(); T sign_bit = pand(should_flipsign, sign_mask); return pxor(sign_bit, x); }