From 8eaa7552fe3147b55d522da910a22fe5aa99d2a5 Mon Sep 17 00:00:00 2001
From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com>
Date: Sun, 8 Mar 2026 14:59:23 -0700
Subject: [PATCH] Fix three flaky tests: packetmath, array_cwise,
 polynomialsolver

libeigen/eigen!2267

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
---
 test/array_cwise.cpp                  | 19 ++++---------------
 test/packetmath.cpp                   | 17 +++++++++++------
 unsupported/test/polynomialsolver.cpp | 22 ++++++++++++++--------
 3 files changed, 29 insertions(+), 29 deletions(-)

diff --git a/test/array_cwise.cpp b/test/array_cwise.cpp
index 64528b928..f73a3bf85 100644
--- a/test/array_cwise.cpp
+++ b/test/array_cwise.cpp
@@ -1006,21 +1006,10 @@ void array_complex(const ArrayType& m) {
 
   std::complex<RealScalar> zero(0.0, 0.0);
   VERIFY((Eigen::isnan)(m1 * zero / zero).all());
-#if EIGEN_COMP_MSVC
-  // msvc complex division is not robust
-  VERIFY((Eigen::isinf)(m4 / RealScalar(0)).all());
-#else
-#if EIGEN_COMP_CLANG
-  // clang's complex division is notoriously broken too
-  if ((numext::isinf)(m4(0, 0) / RealScalar(0))) {
-#endif
-    VERIFY((Eigen::isinf)(m4 / zero).all());
-#if EIGEN_COMP_CLANG
-  } else {
-    VERIFY((Eigen::isinf)(m4.real() / zero.real()).all());
-  }
-#endif
-#endif  // MSVC
+  // Complex division by zero may produce inf or NaN depending on the std::complex
+  // implementation (algebraic formula gives 0/0=NaN for some components). Only require
+  // the result to be non-finite.
+  VERIFY((!(Eigen::isfinite)(m4 / zero)).all());
 
   VERIFY(((Eigen::isfinite)(m1) && (!(Eigen::isfinite)(m1 * zero / zero)) && (!(Eigen::isfinite)(m1 / zero))).all());
 
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index 33fb852a2..b4300b316 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -1221,7 +1221,17 @@ std::enable_if_t<Cond, void> run_ieee_cases(const FunctorT& fun) {
   const Scalar norm_max = (std::numeric_limits<Scalar>::max)();
   const Scalar inf = (std::numeric_limits<Scalar>::infinity)();
   const Scalar nan = (std::numeric_limits<Scalar>::quiet_NaN)();
-  std::vector<Scalar> values{norm_min, Scalar(0), Scalar(1), norm_max, inf, nan};
+  std::vector<Scalar> values{Scalar(0), Scalar(1), norm_max, inf, nan};
+  // On ARM, NEON flush-to-zero mode can flush intermediate subnormal results to zero,
+  // causing functions like sin(norm_min) to return 0 instead of norm_min. Skip norm_min
+  // in that case, along with truly subnormal values.
+  if (!SkipDenorms) {
+    values.push_back(norm_min);
+    if (std::numeric_limits<Scalar>::has_denorm == std::denorm_present) {
+      values.push_back(std::numeric_limits<Scalar>::denorm_min());
+      values.push_back(norm_min / Scalar(2));
+    }
+  }
 
   constexpr int size = PacketSize * 2;
   EIGEN_ALIGN_MAX Scalar data1[size];
@@ -1231,11 +1241,6 @@ std::enable_if_t<Cond, void> run_ieee_cases(const FunctorT& fun) {
     data1[i] = data2[i] = ref[i] = Scalar(0);
   }
 
-  if (Cond && !SkipDenorms && std::numeric_limits<Scalar>::has_denorm == std::denorm_present) {
-    values.push_back(std::numeric_limits<Scalar>::denorm_min());
-    values.push_back(norm_min / Scalar(2));
-  }
-
   for (Scalar abs_value : values) {
     data1[0] = abs_value;
     data1[1] = -data1[0];
diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp
index fc66c4240..1d43f928e 100644
--- a/unsupported/test/polynomialsolver.cpp
+++ b/unsupported/test/polynomialsolver.cpp
@@ -115,7 +115,7 @@ void evalSolverSugarFunction(const POLYNOMIAL& pols, const ROOTS& roots, const R
 
     for (size_t i = 0; i < calc_realRoots.size(); ++i) {
       bool found = false;
-      for (size_t j = 0; j < calc_realRoots.size() && !found; ++j) {
+      for (size_t j = 0; j < real_roots.size() && !found; ++j) {
         if (internal::isApprox(calc_realRoots[i], real_roots[j], psPrec)) {
           found = true;
         }
@@ -178,13 +178,19 @@ void polynomialsolver(int deg) {
   roots_to_monicPolynomial(allRoots, pols);
   evalSolver<Deg_, PolynomialType>(pols);
 
-  cout << "Test sugar" << endl;
-  RealRootsType realRoots = RealRootsType::Random(deg);
-  // sort by ascending absolute value to mitigate precision lost during polynomial expansion
-  std::sort(realRoots.begin(), realRoots.end(),
-            [](RealScalar a, RealScalar b) { return numext::abs(a) < numext::abs(b); });
-  roots_to_monicPolynomial(realRoots, pols);
-  evalSolverSugarFunction<Deg_>(pols, realRoots.template cast<std::complex<RealScalar> >().eval(), realRoots);
+  // The companion matrix eigenvalue approach has limited accuracy for float at
+  // high degrees. The PolynomialSolver documentation itself warns: "With 32bit
+  // (float) floating types this problem shows up frequently." Skip the sugar
+  // function test (which requires exact root matching) for float beyond degree 8.
+  if (deg <= 8 || sizeof(RealScalar) > sizeof(float)) {
+    cout << "Test sugar" << endl;
+    RealRootsType realRoots = RealRootsType::Random(deg);
+    // sort by ascending absolute value to mitigate precision lost during polynomial expansion
+    std::sort(realRoots.begin(), realRoots.end(),
+              [](RealScalar a, RealScalar b) { return numext::abs(a) < numext::abs(b); });
+    roots_to_monicPolynomial(realRoots, pols);
+    evalSolverSugarFunction<Deg_>(pols, realRoots.template cast<std::complex<RealScalar> >().eval(), realRoots);
+  }
 }
 
 EIGEN_DECLARE_TEST(polynomialsolver) {