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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
00026 #define EIGEN_POLYNOMIAL_SOLVER_H
00027 
00041 template< typename _Scalar, int _Deg >
00042 class PolynomialSolverBase
00043 {
00044   public:
00045     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00046 
00047     typedef _Scalar                             Scalar;
00048     typedef typename NumTraits<Scalar>::Real    RealScalar;
00049     typedef std::complex<RealScalar>            RootType;
00050     typedef Matrix<RootType,_Deg,1>             RootsType;
00051 
00052     typedef DenseIndex Index;
00053 
00054   protected:
00055     template< typename OtherPolynomial >
00056     inline void setPolynomial( const OtherPolynomial& poly ){
00057       m_roots.resize(poly.size()); }
00058 
00059   public:
00060     template< typename OtherPolynomial >
00061     inline PolynomialSolverBase( const OtherPolynomial& poly ){
00062       setPolynomial( poly() ); }
00063 
00064     inline PolynomialSolverBase(){}
00065 
00066   public:
00068     inline const RootsType& roots() const { return m_roots; }
00069 
00070   public:
00081     template<typename Stl_back_insertion_sequence>
00082     inline void realRoots( Stl_back_insertion_sequence& bi_seq,
00083         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00084     {
00085       bi_seq.clear();
00086       for(Index i=0; i<m_roots.size(); ++i )
00087       {
00088         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
00089           bi_seq.push_back( m_roots[i].real() ); }
00090       }
00091     }
00092 
00093   protected:
00094     template<typename squaredNormBinaryPredicate>
00095     inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
00096     {
00097       Index res=0;
00098       RealScalar norm2 = internal::abs2( m_roots[0] );
00099       for( Index i=1; i<m_roots.size(); ++i )
00100       {
00101         const RealScalar currNorm2 = internal::abs2( m_roots[i] );
00102         if( pred( currNorm2, norm2 ) ){
00103           res=i; norm2=currNorm2; }
00104       }
00105       return m_roots[res];
00106     }
00107 
00108   public:
00112     inline const RootType& greatestRoot() const
00113     {
00114       std::greater<Scalar> greater;
00115       return selectComplexRoot_withRespectToNorm( greater );
00116     }
00117 
00121     inline const RootType& smallestRoot() const
00122     {
00123       std::less<Scalar> less;
00124       return selectComplexRoot_withRespectToNorm( less );
00125     }
00126 
00127   protected:
00128     template<typename squaredRealPartBinaryPredicate>
00129     inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
00130         squaredRealPartBinaryPredicate& pred,
00131         bool& hasArealRoot,
00132         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00133     {
00134       hasArealRoot = false;
00135       Index res=0;
00136       RealScalar abs2(0);
00137 
00138       for( Index i=0; i<m_roots.size(); ++i )
00139       {
00140         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00141         {
00142           if( !hasArealRoot )
00143           {
00144             hasArealRoot = true;
00145             res = i;
00146             abs2 = m_roots[i].real() * m_roots[i].real();
00147           }
00148           else
00149           {
00150             const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
00151             if( pred( currAbs2, abs2 ) )
00152             {
00153               abs2 = currAbs2;
00154               res = i;
00155             }
00156           }
00157         }
00158         else
00159         {
00160           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00161             res = i; }
00162         }
00163       }
00164       return internal::real_ref(m_roots[res]);
00165     }
00166 
00167 
00168     template<typename RealPartBinaryPredicate>
00169     inline const RealScalar& selectRealRoot_withRespectToRealPart(
00170         RealPartBinaryPredicate& pred,
00171         bool& hasArealRoot,
00172         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00173     {
00174       hasArealRoot = false;
00175       Index res=0;
00176       RealScalar val(0);
00177 
00178       for( Index i=0; i<m_roots.size(); ++i )
00179       {
00180         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00181         {
00182           if( !hasArealRoot )
00183           {
00184             hasArealRoot = true;
00185             res = i;
00186             val = m_roots[i].real();
00187           }
00188           else
00189           {
00190             const RealScalar curr = m_roots[i].real();
00191             if( pred( curr, val ) )
00192             {
00193               val = curr;
00194               res = i;
00195             }
00196           }
00197         }
00198         else
00199         {
00200           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00201             res = i; }
00202         }
00203       }
00204       return internal::real_ref(m_roots[res]);
00205     }
00206 
00207   public:
00222     inline const RealScalar& absGreatestRealRoot(
00223         bool& hasArealRoot,
00224         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00225     {
00226       std::greater<Scalar> greater;
00227       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
00228     }
00229 
00230 
00245     inline const RealScalar& absSmallestRealRoot(
00246         bool& hasArealRoot,
00247         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00248     {
00249       std::less<Scalar> less;
00250       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
00251     }
00252 
00253 
00268     inline const RealScalar& greatestRealRoot(
00269         bool& hasArealRoot,
00270         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00271     {
00272       std::greater<Scalar> greater;
00273       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
00274     }
00275 
00276 
00291     inline const RealScalar& smallestRealRoot(
00292         bool& hasArealRoot,
00293         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00294     {
00295       std::less<Scalar> less;
00296       return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
00297     }
00298 
00299   protected:
00300     RootsType               m_roots;
00301 };
00302 
00303 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
00304   typedef typename BASE::Scalar                 Scalar;       \
00305   typedef typename BASE::RealScalar             RealScalar;   \
00306   typedef typename BASE::RootType               RootType;     \
00307   typedef typename BASE::RootsType              RootsType;
00308 
00309 
00310 
00340 template< typename _Scalar, int _Deg >
00341 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
00342 {
00343   public:
00344     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00345 
00346     typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
00347     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00348 
00349     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
00350     typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
00351 
00352   public:
00354     template< typename OtherPolynomial >
00355     void compute( const OtherPolynomial& poly )
00356     {
00357       assert( Scalar(0) != poly[poly.size()-1] );
00358       internal::companion<Scalar,_Deg> companion( poly );
00359       companion.balance();
00360       m_eigenSolver.compute( companion.denseMatrix() );
00361       m_roots = m_eigenSolver.eigenvalues();
00362     }
00363 
00364   public:
00365     template< typename OtherPolynomial >
00366     inline PolynomialSolver( const OtherPolynomial& poly ){
00367       compute( poly ); }
00368 
00369     inline PolynomialSolver(){}
00370 
00371   protected:
00372     using                   PS_Base::m_roots;
00373     EigenSolverType         m_eigenSolver;
00374 };
00375 
00376 
00377 template< typename _Scalar >
00378 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
00379 {
00380   public:
00381     typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
00382     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00383 
00384   public:
00386     template< typename OtherPolynomial >
00387     void compute( const OtherPolynomial& poly )
00388     {
00389       assert( Scalar(0) != poly[poly.size()-1] );
00390       m_roots[0] = -poly[0]/poly[poly.size()-1];
00391     }
00392 
00393   protected:
00394     using                   PS_Base::m_roots;
00395 };
00396 
00397 #endif // EIGEN_POLYNOMIAL_SOLVER_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:11