1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_POLYNOMIAL_SOLVER_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_POLYNOMIAL_SOLVER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Polynomials_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PolynomialSolverBase. 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Defined to be inherited by polynomial solvers: it provides 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * convenient methods such as 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - real roots, 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - greatest, smallest complex roots, 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - real roots with greatest, smallest absolute real value, 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - greatest, smallest real roots. 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It stores the set of roots as a vector of complexes. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _Scalar, int _Deg > 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PolynomialSolverBase 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::complex<RealScalar> RootType; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<RootType,_Deg,1> RootsType; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef DenseIndex Index; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template< typename OtherPolynomial > 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void setPolynomial( const OtherPolynomial& poly ){ 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_roots.resize(poly.size()); } 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template< typename OtherPolynomial > 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline PolynomialSolverBase( const OtherPolynomial& poly ){ 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath setPolynomial( poly() ); } 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline PolynomialSolverBase(){} 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the complex roots of the polynomial */ 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RootsType& roots() const { return m_roots; } 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Clear and fills the back insertion sequence with the real roots of the polynomial 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * i.e. the real part of the complex roots that have an imaginary part which 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absolute value is smaller than absImaginaryThreshold. 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold takes the dummy_precision associated 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with the _Scalar template parameter of the PolynomialSolver class as the default value. 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] bi_seq : the back insertion sequence (stl concept) 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * number that is considered as real. 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * */ 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Stl_back_insertion_sequence> 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void realRoots( Stl_back_insertion_sequence& bi_seq, 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bi_seq.clear(); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=0; i<m_roots.size(); ++i ) 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){ 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bi_seq.push_back( m_roots[i].real() ); } 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename squaredNormBinaryPredicate> 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index res=0; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar norm2 = internal::abs2( m_roots[0] ); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for( Index i=1; i<m_roots.size(); ++i ) 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar currNorm2 = internal::abs2( m_roots[i] ); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( pred( currNorm2, norm2 ) ){ 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res=i; norm2=currNorm2; } 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_roots[res]; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the complex root with greatest norm. 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RootType& greatestRoot() const 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::greater<Scalar> greater; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectComplexRoot_withRespectToNorm( greater ); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the complex root with smallest norm. 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RootType& smallestRoot() const 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::less<Scalar> less; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectComplexRoot_withRespectToNorm( less ); 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename squaredRealPartBinaryPredicate> 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath squaredRealPartBinaryPredicate& pred, 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hasArealRoot = false; 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index res=0; 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar abs2(0); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for( Index i=0; i<m_roots.size(); ++i ) 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ) 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( !hasArealRoot ) 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hasArealRoot = true; 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath abs2 = m_roots[i].real() * m_roots[i].real(); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( pred( currAbs2, abs2 ) ) 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath abs2 = currAbs2; 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){ 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; } 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::real_ref(m_roots[res]); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename RealPartBinaryPredicate> 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& selectRealRoot_withRespectToRealPart( 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealPartBinaryPredicate& pred, 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hasArealRoot = false; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index res=0; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar val(0); 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for( Index i=0; i<m_roots.size(); ++i ) 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ) 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( !hasArealRoot ) 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hasArealRoot = true; 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath val = m_roots[i].real(); 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar curr = m_roots[i].real(); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( pred( curr, val ) ) 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath val = curr; 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){ 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = i; } 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::real_ref(m_roots[res]); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a real root with greatest absolute magnitude. 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A real root is defined as the real part of a complex root with absolute imaginary 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part smallest than absImaginaryThreshold. 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold takes the dummy_precision associated 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with the _Scalar template parameter of the PolynomialSolver class as the default value. 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If no real root is found the boolean hasArealRoot is set to false and the real part of 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the root with smallest absolute imaginary part is returned instead. 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] hasArealRoot : boolean true if a real root is found according to the 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold criterion, false otherwise. 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * whether or not a root is real. 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& absGreatestRealRoot( 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::greater<Scalar> greater; 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a real root with smallest absolute magnitude. 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A real root is defined as the real part of a complex root with absolute imaginary 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part smallest than absImaginaryThreshold. 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold takes the dummy_precision associated 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with the _Scalar template parameter of the PolynomialSolver class as the default value. 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If no real root is found the boolean hasArealRoot is set to false and the real part of 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the root with smallest absolute imaginary part is returned instead. 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] hasArealRoot : boolean true if a real root is found according to the 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold criterion, false otherwise. 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * whether or not a root is real. 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& absSmallestRealRoot( 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::less<Scalar> less; 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the real root with greatest value. 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A real root is defined as the real part of a complex root with absolute imaginary 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part smallest than absImaginaryThreshold. 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold takes the dummy_precision associated 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with the _Scalar template parameter of the PolynomialSolver class as the default value. 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If no real root is found the boolean hasArealRoot is set to false and the real part of 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the root with smallest absolute imaginary part is returned instead. 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] hasArealRoot : boolean true if a real root is found according to the 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold criterion, false otherwise. 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * whether or not a root is real. 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& greatestRealRoot( 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::greater<Scalar> greater; 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the real root with smallest value. 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * A real root is defined as the real part of a complex root with absolute imaginary 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part smallest than absImaginaryThreshold. 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold takes the dummy_precision associated 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with the _Scalar template parameter of the PolynomialSolver class as the default value. 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If no real root is found the boolean hasArealRoot is set to false and the real part of 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the root with smallest absolute imaginary part is returned instead. 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] hasArealRoot : boolean true if a real root is found according to the 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * absImaginaryThreshold criterion, false otherwise. 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * whether or not a root is real. 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const RealScalar& smallestRealRoot( 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool& hasArealRoot, 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::less<Scalar> less; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RootsType m_roots; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename BASE::Scalar Scalar; \ 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename BASE::RealScalar RealScalar; \ 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename BASE::RootType RootType; \ 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename BASE::RootsType RootsType; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Polynomials_Module 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class PolynomialSolver 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief A polynomial solver 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computes the complex roots of a real polynomial. 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic. 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Notice that the number of polynomial coefficients is _Deg+1. 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class implements a polynomial solver and provides convenient methods such as 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - real roots, 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - greatest, smallest complex roots, 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - real roots with greatest, smallest absolute real value. 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - greatest, smallest real roots. 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules. 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the polynomial to compute its roots. 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This supposes that the complex moduli of the roots are all distinct: e.g. there should 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * be no multiple roots or conjugate roots for instance. 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * With 32bit (float) floating types this problem shows up frequently. 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * However, almost always, correct accuracy is reached even in these cases for 64bit 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (double) floating types and small polynomial degree (<20). 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _Scalar, int _Deg > 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef EigenSolver<CompanionMatrixType> EigenSolverType; 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the complex roots of a new polynomial. */ 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template< typename OtherPolynomial > 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute( const OtherPolynomial& poly ) 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert( Scalar(0) != poly[poly.size()-1] ); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::companion<Scalar,_Deg> companion( poly ); 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath companion.balance(); 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenSolver.compute( companion.denseMatrix() ); 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_roots = m_eigenSolver.eigenvalues(); 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template< typename OtherPolynomial > 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline PolynomialSolver( const OtherPolynomial& poly ){ 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute( poly ); } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline PolynomialSolver(){} 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using PS_Base::m_roots; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolverType m_eigenSolver; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate< typename _Scalar > 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PolynomialSolverBase<_Scalar,1> PS_Base; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Computes the complex roots of a new polynomial. */ 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template< typename OtherPolynomial > 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute( const OtherPolynomial& poly ) 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert( Scalar(0) != poly[poly.size()-1] ); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_roots[0] = -poly[0]/poly[poly.size()-1]; 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using PS_Base::m_roots; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_POLYNOMIAL_SOLVER_H 387