diff options
-rw-r--r-- | unsupported/Eigen/src/Polynomials/Companion.h | 3 | ||||
-rw-r--r-- | unsupported/Eigen/src/Polynomials/PolynomialSolver.h | 30 | ||||
-rw-r--r-- | unsupported/test/polynomialsolver.cpp | 21 |
3 files changed, 46 insertions, 8 deletions
diff --git a/unsupported/Eigen/src/Polynomials/Companion.h b/unsupported/Eigen/src/Polynomials/Companion.h index 126be783b..6ab8f9714 100644 --- a/unsupported/Eigen/src/Polynomials/Companion.h +++ b/unsupported/Eigen/src/Polynomials/Companion.h @@ -75,8 +75,7 @@ class companion void setPolynomial( const VectorType& poly ) { const Index deg = poly.size()-1; - m_monic = Scalar(-1)/poly[deg] * poly.head(deg); - //m_bl_diag.setIdentity( deg-1 ); + m_monic = -poly.head(deg)/poly[deg]; m_bl_diag.setOnes(deg-1); } diff --git a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index 788594247..d4f737134 100644 --- a/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -126,7 +126,7 @@ class PolynomialSolverBase 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 ) { @@ -144,10 +144,10 @@ class PolynomialSolverBase } } } - else + else if(!hasArealRoot) { if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ - res = i; } + res = i;} } } return numext::real_ref(m_roots[res]); @@ -167,7 +167,7 @@ class PolynomialSolverBase 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 ) { @@ -340,6 +340,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexEigenSolver<CompanionMatrixType>, EigenSolver<CompanionMatrixType> >::type EigenSolverType; + typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar> >::type ComplexScalar; public: /** Computes the complex roots of a new polynomial. */ @@ -354,6 +355,27 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> companion.balance(); m_eigenSolver.compute( companion.denseMatrix() ); m_roots = m_eigenSolver.eigenvalues(); + MatrixXcd A = companion.denseMatrix(); + // cleanup noise in imaginary part of real roots: + // if the imaginary part is rather small compared to the real part + // and that cancelling the imaginary part yield a smaller evaluation, + // then it's safe to keep the real part only. + RealScalar coarse_prec = std::pow(4,poly.size()+1)*NumTraits<RealScalar>::epsilon(); + std::cout << coarse_prec << "\n"; + for(Index i = 0; i<m_roots.size(); ++i) + { + if( internal::isMuchSmallerThan(numext::abs(numext::imag(m_roots[i])), + numext::abs(numext::real(m_roots[i])), + coarse_prec) ) + { + ComplexScalar as_real_root = ComplexScalar(numext::real(m_roots[i])); + if( numext::abs(poly_eval(poly, as_real_root)) + <= numext::abs(poly_eval(poly, m_roots[i]))) + { + m_roots[i] = as_real_root; + } + } + } } else if(poly.size () == 2) { diff --git a/unsupported/test/polynomialsolver.cpp b/unsupported/test/polynomialsolver.cpp index 50c74f797..4ff9bda5a 100644 --- a/unsupported/test/polynomialsolver.cpp +++ b/unsupported/test/polynomialsolver.cpp @@ -26,6 +26,16 @@ struct increment_if_fixed_size } } +template<typename PolynomialType> +PolynomialType polyder(const PolynomialType& p) +{ + typedef typename PolynomialType::Scalar Scalar; + PolynomialType res(p.size()); + for(Index i=1; i<p.size(); ++i) + res[i-1] = p[i]*Scalar(i); + res[p.size()-1] = 0.; + return res; +} template<int Deg, typename POLYNOMIAL, typename SOLVER> bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) @@ -44,10 +54,17 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve ) psolve.compute( pols ); const RootsType& roots( psolve.roots() ); EvalRootsType evr( deg ); + POLYNOMIAL pols_der = polyder(pols); + EvalRootsType der( deg ); for( int i=0; i<roots.size(); ++i ){ - evr[i] = std::abs( poly_eval( pols, roots[i] ) ); } + evr[i] = std::abs( poly_eval( pols, roots[i] ) ); + der[i] = numext::maxi(RealScalar(1.), std::abs( poly_eval( pols_der, roots[i] ) )); + } - bool evalToZero = evr.isZero( test_precision<Scalar>() ); + // we need to divide by the magnitude of the derivative because + // with a high derivative is very small error in the value of the root + // yiels a very large error in the polynomial evaluation. + bool evalToZero = (evr.cwiseQuotient(der)).isZero( test_precision<Scalar>() ); if( !evalToZero ) { cerr << "WRONG root: " << endl; |