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       bi_seq.clear();
00073       for(Index i=0; i<m_roots.size(); ++i )
00074       {
00075         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
00076           bi_seq.push_back( m_roots[i].real() ); }
00077       }
00078     }
00079 
00080   protected:
00081     template<typename squaredNormBinaryPredicate>
00082     inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
00083     {
00084       Index res=0;
00085       RealScalar norm2 = internal::abs2( m_roots[0] );
00086       for( Index i=1; i<m_roots.size(); ++i )
00087       {
00088         const RealScalar currNorm2 = internal::abs2( m_roots[i] );
00089         if( pred( currNorm2, norm2 ) ){
00090           res=i; norm2=currNorm2; }
00091       }
00092       return m_roots[res];
00093     }
00094 
00095   public:
00099     inline const RootType& greatestRoot() const
00100     {
00101       std::greater<Scalar> greater;
00102       return selectComplexRoot_withRespectToNorm( greater );
00103     }
00104 
00108     inline const RootType& smallestRoot() const
00109     {
00110       std::less<Scalar> less;
00111       return selectComplexRoot_withRespectToNorm( less );
00112     }
00113 
00114   protected:
00115     template<typename squaredRealPartBinaryPredicate>
00116     inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
00117         squaredRealPartBinaryPredicate& pred,
00118         bool& hasArealRoot,
00119         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00120     {
00121       hasArealRoot = false;
00122       Index res=0;
00123       RealScalar abs2(0);
00124 
00125       for( Index i=0; i<m_roots.size(); ++i )
00126       {
00127         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00128         {
00129           if( !hasArealRoot )
00130           {
00131             hasArealRoot = true;
00132             res = i;
00133             abs2 = m_roots[i].real() * m_roots[i].real();
00134           }
00135           else
00136           {
00137             const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
00138             if( pred( currAbs2, abs2 ) )
00139             {
00140               abs2 = currAbs2;
00141               res = i;
00142             }
00143           }
00144         }
00145         else
00146         {
00147           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00148             res = i; }
00149         }
00150       }
00151       return internal::real_ref(m_roots[res]);
00152     }
00153 
00154 
00155     template<typename RealPartBinaryPredicate>
00156     inline const RealScalar& selectRealRoot_withRespectToRealPart(
00157         RealPartBinaryPredicate& pred,
00158         bool& hasArealRoot,
00159         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00160     {
00161       hasArealRoot = false;
00162       Index res=0;
00163       RealScalar val(0);
00164 
00165       for( Index i=0; i<m_roots.size(); ++i )
00166       {
00167         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00168         {
00169           if( !hasArealRoot )
00170           {
00171             hasArealRoot = true;
00172             res = i;
00173             val = m_roots[i].real();
00174           }
00175           else
00176           {
00177             const RealScalar curr = m_roots[i].real();
00178             if( pred( curr, val ) )
00179             {
00180               val = curr;
00181               res = i;
00182             }
00183           }
00184         }
00185         else
00186         {
00187           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00188             res = i; }
00189         }
00190       }
00191       return internal::real_ref(m_roots[res]);
00192     }
00193 
00194   public:
00209     inline const RealScalar& absGreatestRealRoot(
00210         bool& hasArealRoot,
00211         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00212     {
00213       std::greater<Scalar> greater;
00214       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
00215     }
00216 
00217 
00232     inline const RealScalar& absSmallestRealRoot(
00233         bool& hasArealRoot,
00234         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00235     {
00236       std::less<Scalar> less;
00237       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
00238     }
00239 
00240 
00255     inline const RealScalar& greatestRealRoot(
00256         bool& hasArealRoot,
00257         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00258     {
00259       std::greater<Scalar> greater;
00260       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
00261     }
00262 
00263 
00278     inline const RealScalar& smallestRealRoot(
00279         bool& hasArealRoot,
00280         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00281     {
00282       std::less<Scalar> less;
00283       return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
00284     }
00285 
00286   protected:
00287     RootsType               m_roots;
00288 };
00289 
00290 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
00291   typedef typename BASE::Scalar                 Scalar;       \
00292   typedef typename BASE::RealScalar             RealScalar;   \
00293   typedef typename BASE::RootType               RootType;     \
00294   typedef typename BASE::RootsType              RootsType;
00295 
00296 
00297 
00327 template< typename _Scalar, int _Deg >
00328 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
00329 {
00330   public:
00331     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00332 
00333     typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
00334     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00335 
00336     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
00337     typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
00338 
00339   public:
00341     template< typename OtherPolynomial >
00342     void compute( const OtherPolynomial& poly )
00343     {
00344       assert( Scalar(0) != poly[poly.size()-1] );
00345       internal::companion<Scalar,_Deg> companion( poly );
00346       companion.balance();
00347       m_eigenSolver.compute( companion.denseMatrix() );
00348       m_roots = m_eigenSolver.eigenvalues();
00349     }
00350 
00351   public:
00352     template< typename OtherPolynomial >
00353     inline PolynomialSolver( const OtherPolynomial& poly ){
00354       compute( poly ); }
00355 
00356     inline PolynomialSolver(){}
00357 
00358   protected:
00359     using                   PS_Base::m_roots;
00360     EigenSolverType         m_eigenSolver;
00361 };
00362 
00363 
00364 template< typename _Scalar >
00365 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
00366 {
00367   public:
00368     typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
00369     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00370 
00371   public:
00373     template< typename OtherPolynomial >
00374     void compute( const OtherPolynomial& poly )
00375     {
00376       assert( Scalar(0) != poly[poly.size()-1] );
00377       m_roots[0] = -poly[0]/poly[poly.size()-1];
00378     }
00379 
00380   protected:
00381     using                   PS_Base::m_roots;
00382 };
00383 
00384 } // end namespace Eigen
00385 
00386 #endif // EIGEN_POLYNOMIAL_SOLVER_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:44