blob: 64cab12045cf9c06e8107af70b79658416ed54aa [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2010 Manuel Yguel // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_POLYNOMIAL_SOLVER_H #define EIGEN_POLYNOMIAL_SOLVER_H #include "./InternalHeaderCheck.h" namespace Eigen { /** \ingroup Polynomials_Module * \class PolynomialSolverBase. * * \brief Defined to be inherited by polynomial solvers: it provides * convenient methods such as * - real roots, * - greatest, smallest complex roots, * - real roots with greatest, smallest absolute real value, * - greatest, smallest real roots. * * It stores the set of roots as a vector of complexes. * */ template< typename Scalar_, int Deg_ > class PolynomialSolverBase { public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar_,Deg_==Dynamic ? Dynamic : Deg_) typedef Scalar_ Scalar; typedef typename NumTraits::Real RealScalar; typedef std::complex RootType; typedef Matrix RootsType; typedef DenseIndex Index; protected: template< typename OtherPolynomial > inline void setPolynomial( const OtherPolynomial& poly ){ m_roots.resize(poly.size()-1); } public: template< typename OtherPolynomial > inline PolynomialSolverBase( const OtherPolynomial& poly ){ setPolynomial( poly() ); } inline PolynomialSolverBase(){} public: /** \returns the complex roots of the polynomial */ inline const RootsType& roots() const { return m_roots; } public: /** Clear and fills the back insertion sequence with the real roots of the polynomial * i.e. the real part of the complex roots that have an imaginary part which * absolute value is smaller than absImaginaryThreshold. * absImaginaryThreshold takes the dummy_precision associated * with the Scalar_ template parameter of the PolynomialSolver class as the default value. * * \param[out] bi_seq : the back insertion sequence (stl concept) * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex * number that is considered as real. * */ template inline void realRoots( Stl_back_insertion_sequence& bi_seq, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { using std::abs; bi_seq.clear(); for(Index i=0; i inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const { Index res=0; RealScalar norm2 = numext::abs2( m_roots[0] ); for( Index i=1; i greater; return selectComplexRoot_withRespectToNorm( greater ); } /** * \returns the complex root with smallest norm. */ inline const RootType& smallestRoot() const { std::less less; return selectComplexRoot_withRespectToNorm( less ); } protected: template inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( squaredRealPartBinaryPredicate& pred, bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { using std::abs; hasArealRoot = false; Index res=0; RealScalar abs2(0); for( Index i=0; i inline const RealScalar& selectRealRoot_withRespectToRealPart( RealPartBinaryPredicate& pred, bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { using std::abs; hasArealRoot = false; Index res=0; RealScalar val(0); for( Index i=0; i::dummy_precision() ) const { std::greater greater; return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); } /** * \returns a real root with smallest absolute magnitude. * A real root is defined as the real part of a complex root with absolute imaginary * part smallest than absImaginaryThreshold. * absImaginaryThreshold takes the dummy_precision associated * with the Scalar_ template parameter of the PolynomialSolver class as the default value. * If no real root is found the boolean hasArealRoot is set to false and the real part of * the root with smallest absolute imaginary part is returned instead. * * \param[out] hasArealRoot : boolean true if a real root is found according to the * absImaginaryThreshold criterion, false otherwise. * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide * whether or not a root is real. */ inline const RealScalar& absSmallestRealRoot( bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { std::less less; return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); } /** * \returns the real root with greatest value. * A real root is defined as the real part of a complex root with absolute imaginary * part smallest than absImaginaryThreshold. * absImaginaryThreshold takes the dummy_precision associated * with the Scalar_ template parameter of the PolynomialSolver class as the default value. * If no real root is found the boolean hasArealRoot is set to false and the real part of * the root with smallest absolute imaginary part is returned instead. * * \param[out] hasArealRoot : boolean true if a real root is found according to the * absImaginaryThreshold criterion, false otherwise. * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide * whether or not a root is real. */ inline const RealScalar& greatestRealRoot( bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { std::greater greater; return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); } /** * \returns the real root with smallest value. * A real root is defined as the real part of a complex root with absolute imaginary * part smallest than absImaginaryThreshold. * absImaginaryThreshold takes the dummy_precision associated * with the Scalar_ template parameter of the PolynomialSolver class as the default value. * If no real root is found the boolean hasArealRoot is set to false and the real part of * the root with smallest absolute imaginary part is returned instead. * * \param[out] hasArealRoot : boolean true if a real root is found according to the * absImaginaryThreshold criterion, false otherwise. * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide * whether or not a root is real. */ inline const RealScalar& smallestRealRoot( bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits::dummy_precision() ) const { std::less less; return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); } protected: RootsType m_roots; }; #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ typedef typename BASE::Scalar Scalar; \ typedef typename BASE::RealScalar RealScalar; \ typedef typename BASE::RootType RootType; \ typedef typename BASE::RootsType RootsType; /** \ingroup Polynomials_Module * * \class PolynomialSolver * * \brief A polynomial solver * * Computes the complex roots of a real polynomial. * * \param Scalar_ the scalar type, i.e., the type of the polynomial coefficients * \param Deg_ the degree of the polynomial, can be a compile time value or Dynamic. * Notice that the number of polynomial coefficients is Deg_+1. * * This class implements a polynomial solver and provides convenient methods such as * - real roots, * - greatest, smallest complex roots, * - real roots with greatest, smallest absolute real value. * - greatest, smallest real roots. * * WARNING: this polynomial solver is experimental, part of the unsupported Eigen modules. * * * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of * the polynomial to compute its roots. * This supposes that the complex moduli of the roots are all distinct: e.g. there should * be no multiple roots or conjugate roots for instance. * With 32bit (float) floating types this problem shows up frequently. * However, almost always, correct accuracy is reached even in these cases for 64bit * (double) floating types and small polynomial degree (<20). */ template class PolynomialSolver : public PolynomialSolverBase { public: EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar_,Deg_==Dynamic ? Dynamic : Deg_) typedef PolynomialSolverBase PS_Base; EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) typedef Matrix CompanionMatrixType; typedef std::conditional_t::IsComplex, ComplexEigenSolver, EigenSolver > EigenSolverType; typedef std::conditional_t::IsComplex, Scalar, std::complex > ComplexScalar; public: /** Computes the complex roots of a new polynomial. */ template< typename OtherPolynomial > void compute( const OtherPolynomial& poly ) { eigen_assert( Scalar(0) != poly[poly.size()-1] ); eigen_assert( poly.size() > 1 ); if(poly.size() > 2 ) { internal::companion companion( poly ); companion.balance(); m_eigenSolver.compute( companion.denseMatrix() ); m_roots = m_eigenSolver.eigenvalues(); // 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 = RealScalar(std::pow(4,poly.size()+1))*NumTraits::epsilon(); for(Index i = 0; i inline PolynomialSolver( const OtherPolynomial& poly ){ compute( poly ); } inline PolynomialSolver(){} protected: using PS_Base::m_roots; EigenSolverType m_eigenSolver; }; template< typename Scalar_ > class PolynomialSolver : public PolynomialSolverBase { public: typedef PolynomialSolverBase PS_Base; EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) public: /** Computes the complex roots of a new polynomial. */ template< typename OtherPolynomial > void compute( const OtherPolynomial& poly ) { eigen_assert( poly.size() == 2 ); eigen_assert( Scalar(0) != poly[1] ); m_roots[0] = -poly[0]/poly[1]; } public: template< typename OtherPolynomial > inline PolynomialSolver( const OtherPolynomial& poly ){ compute( poly ); } inline PolynomialSolver(){} protected: using PS_Base::m_roots; }; } // end namespace Eigen #endif // EIGEN_POLYNOMIAL_SOLVER_H