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