PolynomialSolver.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
00011 #define EIGEN_POLYNOMIAL_SOLVER_H
00012 
00013 namespace Eigen { 
00014 
00028 template< typename _Scalar, int _Deg >
00029 class PolynomialSolverBase
00030 {
00031   public:
00032     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00033 
00034     typedef _Scalar                             Scalar;
00035     typedef typename NumTraits<Scalar>::Real    RealScalar;
00036     typedef std::complex<RealScalar>            RootType;
00037     typedef Matrix<RootType,_Deg,1>             RootsType;
00038 
00039     typedef DenseIndex Index;
00040 
00041   protected:
00042     template< typename OtherPolynomial >
00043     inline void setPolynomial( const OtherPolynomial& poly ){
00044       m_roots.resize(poly.size()); }
00045 
00046   public:
00047     template< typename OtherPolynomial >
00048     inline PolynomialSolverBase( const OtherPolynomial& poly ){
00049       setPolynomial( poly() ); }
00050 
00051     inline PolynomialSolverBase(){}
00052 
00053   public:
00055     inline const RootsType& roots() const { return m_roots; }
00056 
00057   public:
00068     template<typename Stl_back_insertion_sequence>
00069     inline void realRoots( Stl_back_insertion_sequence& bi_seq,
00070         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00071     {
00072       using std::abs;
00073       bi_seq.clear();
00074       for(Index i=0; i<m_roots.size(); ++i )
00075       {
00076         if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){
00077           bi_seq.push_back( m_roots[i].real() ); }
00078       }
00079     }
00080 
00081   protected:
00082     template<typename squaredNormBinaryPredicate>
00083     inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
00084     {
00085       Index res=0;
00086       RealScalar norm2 = numext::abs2( m_roots[0] );
00087       for( Index i=1; i<m_roots.size(); ++i )
00088       {
00089         const RealScalar currNorm2 = numext::abs2( m_roots[i] );
00090         if( pred( currNorm2, norm2 ) ){
00091           res=i; norm2=currNorm2; }
00092       }
00093       return m_roots[res];
00094     }
00095 
00096   public:
00100     inline const RootType& greatestRoot() const
00101     {
00102       std::greater<Scalar> greater;
00103       return selectComplexRoot_withRespectToNorm( greater );
00104     }
00105 
00109     inline const RootType& smallestRoot() const
00110     {
00111       std::less<Scalar> less;
00112       return selectComplexRoot_withRespectToNorm( less );
00113     }
00114 
00115   protected:
00116     template<typename squaredRealPartBinaryPredicate>
00117     inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
00118         squaredRealPartBinaryPredicate& pred,
00119         bool& hasArealRoot,
00120         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00121     {
00122       using std::abs;
00123       hasArealRoot = false;
00124       Index res=0;
00125       RealScalar abs2(0);
00126 
00127       for( Index i=0; i<m_roots.size(); ++i )
00128       {
00129         if( abs( m_roots[i].imag() ) < absImaginaryThreshold )
00130         {
00131           if( !hasArealRoot )
00132           {
00133             hasArealRoot = true;
00134             res = i;
00135             abs2 = m_roots[i].real() * m_roots[i].real();
00136           }
00137           else
00138           {
00139             const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
00140             if( pred( currAbs2, abs2 ) )
00141             {
00142               abs2 = currAbs2;
00143               res = i;
00144             }
00145           }
00146         }
00147         else
00148         {
00149           if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
00150             res = i; }
00151         }
00152       }
00153       return numext::real_ref(m_roots[res]);
00154     }
00155 
00156 
00157     template<typename RealPartBinaryPredicate>
00158     inline const RealScalar& selectRealRoot_withRespectToRealPart(
00159         RealPartBinaryPredicate& pred,
00160         bool& hasArealRoot,
00161         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00162     {
00163       using std::abs;
00164       hasArealRoot = false;
00165       Index res=0;
00166       RealScalar val(0);
00167 
00168       for( Index i=0; i<m_roots.size(); ++i )
00169       {
00170         if( abs( m_roots[i].imag() ) < absImaginaryThreshold )
00171         {
00172           if( !hasArealRoot )
00173           {
00174             hasArealRoot = true;
00175             res = i;
00176             val = m_roots[i].real();
00177           }
00178           else
00179           {
00180             const RealScalar curr = m_roots[i].real();
00181             if( pred( curr, val ) )
00182             {
00183               val = curr;
00184               res = i;
00185             }
00186           }
00187         }
00188         else
00189         {
00190           if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
00191             res = i; }
00192         }
00193       }
00194       return numext::real_ref(m_roots[res]);
00195     }
00196 
00197   public:
00212     inline const RealScalar& absGreatestRealRoot(
00213         bool& hasArealRoot,
00214         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00215     {
00216       std::greater<Scalar> greater;
00217       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
00218     }
00219 
00220 
00235     inline const RealScalar& absSmallestRealRoot(
00236         bool& hasArealRoot,
00237         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00238     {
00239       std::less<Scalar> less;
00240       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
00241     }
00242 
00243 
00258     inline const RealScalar& greatestRealRoot(
00259         bool& hasArealRoot,
00260         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00261     {
00262       std::greater<Scalar> greater;
00263       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
00264     }
00265 
00266 
00281     inline const RealScalar& smallestRealRoot(
00282         bool& hasArealRoot,
00283         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00284     {
00285       std::less<Scalar> less;
00286       return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
00287     }
00288 
00289   protected:
00290     RootsType               m_roots;
00291 };
00292 
00293 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
00294   typedef typename BASE::Scalar                 Scalar;       \
00295   typedef typename BASE::RealScalar             RealScalar;   \
00296   typedef typename BASE::RootType               RootType;     \
00297   typedef typename BASE::RootsType              RootsType;
00298 
00299 
00300 
00330 template< typename _Scalar, int _Deg >
00331 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
00332 {
00333   public:
00334     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00335 
00336     typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
00337     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00338 
00339     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
00340     typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
00341 
00342   public:
00344     template< typename OtherPolynomial >
00345     void compute( const OtherPolynomial& poly )
00346     {
00347       eigen_assert( Scalar(0) != poly[poly.size()-1] );
00348       internal::companion<Scalar,_Deg> companion( poly );
00349       companion.balance();
00350       m_eigenSolver.compute( companion.denseMatrix() );
00351       m_roots = m_eigenSolver.eigenvalues();
00352     }
00353 
00354   public:
00355     template< typename OtherPolynomial >
00356     inline PolynomialSolver( const OtherPolynomial& poly ){
00357       compute( poly ); }
00358 
00359     inline PolynomialSolver(){}
00360 
00361   protected:
00362     using                   PS_Base::m_roots;
00363     EigenSolverType         m_eigenSolver;
00364 };
00365 
00366 
00367 template< typename _Scalar >
00368 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
00369 {
00370   public:
00371     typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
00372     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00373 
00374   public:
00376     template< typename OtherPolynomial >
00377     void compute( const OtherPolynomial& poly )
00378     {
00379       eigen_assert( Scalar(0) != poly[poly.size()-1] );
00380       m_roots[0] = -poly[0]/poly[poly.size()-1];
00381     }
00382 
00383   protected:
00384     using                   PS_Base::m_roots;
00385 };
00386 
00387 } // end namespace Eigen
00388 
00389 #endif // EIGEN_POLYNOMIAL_SOLVER_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:42