PolynomialSolver.h
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 
00028 namespace Eigen { 
00029 
00043 template< typename _Scalar, int _Deg >
00044 class PolynomialSolverBase
00045 {
00046   public:
00047     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00048 
00049     typedef _Scalar                             Scalar;
00050     typedef typename NumTraits<Scalar>::Real    RealScalar;
00051     typedef std::complex<RealScalar>            RootType;
00052     typedef Matrix<RootType,_Deg,1>             RootsType;
00053 
00054     typedef DenseIndex Index;
00055 
00056   protected:
00057     template< typename OtherPolynomial >
00058     inline void setPolynomial( const OtherPolynomial& poly ){
00059       m_roots.resize(poly.size()); }
00060 
00061   public:
00062     template< typename OtherPolynomial >
00063     inline PolynomialSolverBase( const OtherPolynomial& poly ){
00064       setPolynomial( poly() ); }
00065 
00066     inline PolynomialSolverBase(){}
00067 
00068   public:
00070     inline const RootsType& roots() const { return m_roots; }
00071 
00072   public:
00083     template<typename Stl_back_insertion_sequence>
00084     inline void realRoots( Stl_back_insertion_sequence& bi_seq,
00085         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00086     {
00087       bi_seq.clear();
00088       for(Index i=0; i<m_roots.size(); ++i )
00089       {
00090         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
00091           bi_seq.push_back( m_roots[i].real() ); }
00092       }
00093     }
00094 
00095   protected:
00096     template<typename squaredNormBinaryPredicate>
00097     inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
00098     {
00099       Index res=0;
00100       RealScalar norm2 = internal::abs2( m_roots[0] );
00101       for( Index i=1; i<m_roots.size(); ++i )
00102       {
00103         const RealScalar currNorm2 = internal::abs2( m_roots[i] );
00104         if( pred( currNorm2, norm2 ) ){
00105           res=i; norm2=currNorm2; }
00106       }
00107       return m_roots[res];
00108     }
00109 
00110   public:
00114     inline const RootType& greatestRoot() const
00115     {
00116       std::greater<Scalar> greater;
00117       return selectComplexRoot_withRespectToNorm( greater );
00118     }
00119 
00123     inline const RootType& smallestRoot() const
00124     {
00125       std::less<Scalar> less;
00126       return selectComplexRoot_withRespectToNorm( less );
00127     }
00128 
00129   protected:
00130     template<typename squaredRealPartBinaryPredicate>
00131     inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
00132         squaredRealPartBinaryPredicate& pred,
00133         bool& hasArealRoot,
00134         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00135     {
00136       hasArealRoot = false;
00137       Index res=0;
00138       RealScalar abs2(0);
00139 
00140       for( Index i=0; i<m_roots.size(); ++i )
00141       {
00142         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00143         {
00144           if( !hasArealRoot )
00145           {
00146             hasArealRoot = true;
00147             res = i;
00148             abs2 = m_roots[i].real() * m_roots[i].real();
00149           }
00150           else
00151           {
00152             const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
00153             if( pred( currAbs2, abs2 ) )
00154             {
00155               abs2 = currAbs2;
00156               res = i;
00157             }
00158           }
00159         }
00160         else
00161         {
00162           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00163             res = i; }
00164         }
00165       }
00166       return internal::real_ref(m_roots[res]);
00167     }
00168 
00169 
00170     template<typename RealPartBinaryPredicate>
00171     inline const RealScalar& selectRealRoot_withRespectToRealPart(
00172         RealPartBinaryPredicate& pred,
00173         bool& hasArealRoot,
00174         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00175     {
00176       hasArealRoot = false;
00177       Index res=0;
00178       RealScalar val(0);
00179 
00180       for( Index i=0; i<m_roots.size(); ++i )
00181       {
00182         if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
00183         {
00184           if( !hasArealRoot )
00185           {
00186             hasArealRoot = true;
00187             res = i;
00188             val = m_roots[i].real();
00189           }
00190           else
00191           {
00192             const RealScalar curr = m_roots[i].real();
00193             if( pred( curr, val ) )
00194             {
00195               val = curr;
00196               res = i;
00197             }
00198           }
00199         }
00200         else
00201         {
00202           if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
00203             res = i; }
00204         }
00205       }
00206       return internal::real_ref(m_roots[res]);
00207     }
00208 
00209   public:
00224     inline const RealScalar& absGreatestRealRoot(
00225         bool& hasArealRoot,
00226         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00227     {
00228       std::greater<Scalar> greater;
00229       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
00230     }
00231 
00232 
00247     inline const RealScalar& absSmallestRealRoot(
00248         bool& hasArealRoot,
00249         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00250     {
00251       std::less<Scalar> less;
00252       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
00253     }
00254 
00255 
00270     inline const RealScalar& greatestRealRoot(
00271         bool& hasArealRoot,
00272         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00273     {
00274       std::greater<Scalar> greater;
00275       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
00276     }
00277 
00278 
00293     inline const RealScalar& smallestRealRoot(
00294         bool& hasArealRoot,
00295         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
00296     {
00297       std::less<Scalar> less;
00298       return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
00299     }
00300 
00301   protected:
00302     RootsType               m_roots;
00303 };
00304 
00305 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
00306   typedef typename BASE::Scalar                 Scalar;       \
00307   typedef typename BASE::RealScalar             RealScalar;   \
00308   typedef typename BASE::RootType               RootType;     \
00309   typedef typename BASE::RootsType              RootsType;
00310 
00311 
00312 
00342 template< typename _Scalar, int _Deg >
00343 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
00344 {
00345   public:
00346     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00347 
00348     typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
00349     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00350 
00351     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
00352     typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
00353 
00354   public:
00356     template< typename OtherPolynomial >
00357     void compute( const OtherPolynomial& poly )
00358     {
00359       assert( Scalar(0) != poly[poly.size()-1] );
00360       internal::companion<Scalar,_Deg> companion( poly );
00361       companion.balance();
00362       m_eigenSolver.compute( companion.denseMatrix() );
00363       m_roots = m_eigenSolver.eigenvalues();
00364     }
00365 
00366   public:
00367     template< typename OtherPolynomial >
00368     inline PolynomialSolver( const OtherPolynomial& poly ){
00369       compute( poly ); }
00370 
00371     inline PolynomialSolver(){}
00372 
00373   protected:
00374     using                   PS_Base::m_roots;
00375     EigenSolverType         m_eigenSolver;
00376 };
00377 
00378 
00379 template< typename _Scalar >
00380 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
00381 {
00382   public:
00383     typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
00384     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
00385 
00386   public:
00388     template< typename OtherPolynomial >
00389     void compute( const OtherPolynomial& poly )
00390     {
00391       assert( Scalar(0) != poly[poly.size()-1] );
00392       m_roots[0] = -poly[0]/poly[poly.size()-1];
00393     }
00394 
00395   protected:
00396     using                   PS_Base::m_roots;
00397 };
00398 
00399 } // end namespace Eigen
00400 
00401 #endif // EIGEN_POLYNOMIAL_SOLVER_H