Mirrors/eigen
0
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-01-18 17:31:19 +01:00
Code Packages Releases Activity

fix polynomialsolver test failures

Browse Source 
libeigen/eigen!2104
This commit is contained in:
Charles Schlosser
2026-01-05 05:19:49 +00:00
parent 711118b747
commit d90a0534be
2 changed files with 18 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
unsupported/Eigen/src/Polynomials/PolynomialSolver.h
View File

@@ -82,8 +82,8 @@ class PolynomialSolverBase {
} }
protected: protected:
template <typename squaredNormBinaryPredicate> template <typename Predicate>
inline const RootType& selectComplexRoot_withRespectToNorm(squaredNormBinaryPredicate& pred) const { inline const RootType& selectComplexRoot_withRespectToNorm(Predicate& pred) const {
Index res = 0; Index res = 0;
RealScalar norm2 = numext::abs2(m_roots[0]); RealScalar norm2 = numext::abs2(m_roots[0]);
for (Index i = 1; i < m_roots.size(); ++i) { for (Index i = 1; i < m_roots.size(); ++i) {
@@ -114,25 +114,25 @@ class PolynomialSolverBase {
} }
protected: protected:
template <typename squaredRealPartBinaryPredicate> template <typename Predicate>
inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
squaredRealPartBinaryPredicate& pred, bool& hasArealRoot, Predicate& pred, bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision()) const { const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision()) const {
using std::abs; using std::abs;
hasArealRoot = false; hasArealRoot = false;
Index res = 0; Index res = 0;
RealScalar abs2(0); RealScalar val(0);
for (Index i = 0; i < m_roots.size(); ++i) { for (Index i = 0; i < m_roots.size(); ++i) {
if (abs(m_roots[i].imag()) <= absImaginaryThreshold) { if (abs(m_roots[i].imag()) <= absImaginaryThreshold) {
if (!hasArealRoot) { if (!hasArealRoot) {
hasArealRoot = true; hasArealRoot = true;
res = i; res = i;
abs2 = m_roots[i].real() * m_roots[i].real(); val = abs(m_roots[i].real());
} else { } else {
const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); const RealScalar curr = abs(m_roots[i].real());
if (pred(currAbs2, abs2)) { if (pred(curr, val)) {
abs2 = currAbs2; val = curr;
res = i; res = i;
} }
} }
@@ -145,9 +145,9 @@ class PolynomialSolverBase {
return numext::real_ref(m_roots[res]); return numext::real_ref(m_roots[res]);
} }
template <typename RealPartBinaryPredicate> template <typename Predicate>
inline const RealScalar& selectRealRoot_withRespectToRealPart( inline const RealScalar& selectRealRoot_withRespectToRealPart(
RealPartBinaryPredicate& pred, bool& hasArealRoot, Predicate& pred, bool& hasArealRoot,
const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision()) const { const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision()) const {
using std::abs; using std::abs;
hasArealRoot = false; hasArealRoot = false;

11
unsupported/test/polynomialsolver.cpp
View File

@@ -131,28 +131,28 @@ void evalSolverSugarFunction(const POLYNOMIAL& pols, const ROOTS& roots, const R
bool hasRealRoot; bool hasRealRoot;
// Test absGreatestRealRoot // Test absGreatestRealRoot
RealScalar r = psolve.absGreatestRealRoot(hasRealRoot); RealScalar r = psolve.absGreatestRealRoot(hasRealRoot, test_precision<RealScalar>());
VERIFY(hasRealRoot == (real_roots.size() > 0)); VERIFY(hasRealRoot == (real_roots.size() > 0));
if (hasRealRoot) { if (hasRealRoot) {
VERIFY(internal::isApprox(real_roots.array().abs().maxCoeff(), abs(r), psPrec)); VERIFY(internal::isApprox(real_roots.array().abs().maxCoeff(), abs(r), psPrec));
} }
// Test absSmallestRealRoot // Test absSmallestRealRoot
r = psolve.absSmallestRealRoot(hasRealRoot); r = psolve.absSmallestRealRoot(hasRealRoot, test_precision<RealScalar>());
VERIFY(hasRealRoot == (real_roots.size() > 0)); VERIFY(hasRealRoot == (real_roots.size() > 0));
if (hasRealRoot) { if (hasRealRoot) {
VERIFY(internal::isApprox(real_roots.array().abs().minCoeff(), abs(r), psPrec)); VERIFY(internal::isApprox(real_roots.array().abs().minCoeff(), abs(r), psPrec));
} }
// Test greatestRealRoot // Test greatestRealRoot
r = psolve.greatestRealRoot(hasRealRoot); r = psolve.greatestRealRoot(hasRealRoot, test_precision<RealScalar>());
VERIFY(hasRealRoot == (real_roots.size() > 0)); VERIFY(hasRealRoot == (real_roots.size() > 0));
if (hasRealRoot) { if (hasRealRoot) {
VERIFY(internal::isApprox(real_roots.array().maxCoeff(), r, psPrec)); VERIFY(internal::isApprox(real_roots.array().maxCoeff(), r, psPrec));
} }
// Test smallestRealRoot // Test smallestRealRoot
r = psolve.smallestRealRoot(hasRealRoot); r = psolve.smallestRealRoot(hasRealRoot, test_precision<RealScalar>());
VERIFY(hasRealRoot == (real_roots.size() > 0)); VERIFY(hasRealRoot == (real_roots.size() > 0));
if (hasRealRoot) { if (hasRealRoot) {
VERIFY(internal::isApprox(real_roots.array().minCoeff(), r, psPrec)); VERIFY(internal::isApprox(real_roots.array().minCoeff(), r, psPrec));
@@ -180,6 +180,9 @@ void polynomialsolver(int deg) {
cout << "Test sugar" << endl; cout << "Test sugar" << endl;
RealRootsType realRoots = RealRootsType::Random(deg); 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); roots_to_monicPolynomial(realRoots, pols);
evalSolverSugarFunction<Deg_>(pols, realRoots.template cast<std::complex<RealScalar> >().eval(), realRoots); evalSolverSugarFunction<Deg_>(pols, realRoots.template cast<std::complex<RealScalar> >().eval(), realRoots);
} }