| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> |
| // |
| // 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 |
| |
| 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<Scalar>::Real RealScalar; |
| typedef std::complex<RealScalar> RootType; |
| typedef Matrix<RootType,_Deg,1> 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<typename Stl_back_insertion_sequence> |
| inline void realRoots( Stl_back_insertion_sequence& bi_seq, |
| const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const |
| { |
| using std::abs; |
| bi_seq.clear(); |
| for(Index i=0; i<m_roots.size(); ++i ) |
| { |
| if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){ |
| bi_seq.push_back( m_roots[i].real() ); } |
| } |
| } |
| |
| protected: |
| template<typename squaredNormBinaryPredicate> |
| inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const |
| { |
| Index res=0; |
| RealScalar norm2 = numext::abs2( m_roots[0] ); |
| for( Index i=1; i<m_roots.size(); ++i ) |
| { |
| const RealScalar currNorm2 = numext::abs2( m_roots[i] ); |
| if( pred( currNorm2, norm2 ) ){ |
| res=i; norm2=currNorm2; } |
| } |
| return m_roots[res]; |
| } |
| |
| public: |
| /** |
| * \returns the complex root with greatest norm. |
| */ |
| inline const RootType& greatestRoot() const |
| { |
| std::greater<RealScalar> greater; |
| return selectComplexRoot_withRespectToNorm( greater ); |
| } |
| |
| /** |
| * \returns the complex root with smallest norm. |
| */ |
| inline const RootType& smallestRoot() const |
| { |
| std::less<RealScalar> less; |
| return selectComplexRoot_withRespectToNorm( less ); |
| } |
| |
| protected: |
| template<typename squaredRealPartBinaryPredicate> |
| inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( |
| squaredRealPartBinaryPredicate& pred, |
| bool& hasArealRoot, |
| const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const |
| { |
| using std::abs; |
| hasArealRoot = false; |
| Index res=0; |
| RealScalar abs2(0); |
| |
| for( Index i=0; i<m_roots.size(); ++i ) |
| { |
| if( abs( m_roots[i].imag() ) <= absImaginaryThreshold ) |
| { |
| if( !hasArealRoot ) |
| { |
| hasArealRoot = true; |
| res = i; |
| abs2 = m_roots[i].real() * m_roots[i].real(); |
| } |
| else |
| { |
| const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); |
| if( pred( currAbs2, abs2 ) ) |
| { |
| abs2 = currAbs2; |
| res = i; |
| } |
| } |
| } |
| else if(!hasArealRoot) |
| { |
| if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ |
| res = i;} |
| } |
| } |
| return numext::real_ref(m_roots[res]); |
| } |
| |
| |
| template<typename RealPartBinaryPredicate> |
| inline const RealScalar& selectRealRoot_withRespectToRealPart( |
| RealPartBinaryPredicate& pred, |
| bool& hasArealRoot, |
| const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const |
| { |
| using std::abs; |
| hasArealRoot = false; |
| Index res=0; |
| RealScalar val(0); |
| |
| for( Index i=0; i<m_roots.size(); ++i ) |
| { |
| if( abs( m_roots[i].imag() ) <= absImaginaryThreshold ) |
| { |
| if( !hasArealRoot ) |
| { |
| hasArealRoot = true; |
| res = i; |
| val = m_roots[i].real(); |
| } |
| else |
| { |
| const RealScalar curr = m_roots[i].real(); |
| if( pred( curr, val ) ) |
| { |
| val = curr; |
| res = i; |
| } |
| } |
| } |
| else |
| { |
| if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ |
| res = i; } |
| } |
| } |
| return numext::real_ref(m_roots[res]); |
| } |
| |
| public: |
| /** |
| * \returns a real root with greatest 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& absGreatestRealRoot( |
| bool& hasArealRoot, |
| const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const |
| { |
| std::greater<RealScalar> 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<Scalar>::dummy_precision() ) const |
| { |
| std::less<RealScalar> 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<Scalar>::dummy_precision() ) const |
| { |
| std::greater<RealScalar> 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<Scalar>::dummy_precision() ) const |
| { |
| std::less<RealScalar> 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<typename _Scalar, int _Deg> |
| class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> |
| { |
| public: |
| EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) |
| |
| typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; |
| EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) |
| |
| typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; |
| 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. */ |
| 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<Scalar,_Deg> 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<RealScalar>::epsilon(); |
| 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) |
| { |
| m_roots.resize(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; |
| EigenSolverType m_eigenSolver; |
| }; |
| |
| |
| template< typename _Scalar > |
| class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> |
| { |
| public: |
| typedef PolynomialSolverBase<_Scalar,1> 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 |
