EIC code displayed by LXR master/​include/​eigen3/​unsupported/​Eigen/​src/​Polynomials/​PolynomialSolver.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 09:57:06

 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
 #define EIGEN_POLYNOMIAL_SOLVER_H
 
 namespace Eigen { 
 
 /** \ingroup Polynomials_Module
  *  \class PolynomialSolverBase.
  *
  * \brief Defined to be inherited by polynomial solvers: it provides
  * convenient methods such as
  *  - real roots,
  *  - greatest, smallest complex roots,
  *  - real roots with greatest, smallest absolute real value,
  *  - greatest, smallest real roots.
  *
  * It stores the set of roots as a vector of complexes.
  *
  */
 template< typename _Scalar, int _Deg >
 class PolynomialSolverBase
 {
   public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
 
     typedef _Scalar                             Scalar;
     typedef typename NumTraits<Scalar>::Real    RealScalar;
     typedef std::complex<RealScalar>            RootType;
     typedef Matrix<RootType,_Deg,1>             RootsType;
 
     typedef DenseIndex Index;
 
   protected:
     template< typename OtherPolynomial >
     inline void setPolynomial( const OtherPolynomial& poly ){
       m_roots.resize(poly.size()-1); }
 
   public:
     template< typename OtherPolynomial >
     inline PolynomialSolverBase( const OtherPolynomial& poly ){
       setPolynomial( poly() ); }
 
     inline PolynomialSolverBase(){}
 
   public:
     /** \returns the complex roots of the polynomial */
     inline const RootsType& roots() const { return m_roots; }
 
   public:
     /** Clear and fills the back insertion sequence with the real roots of the polynomial
      * i.e. the real part of the complex roots that have an imaginary part which
      * absolute value is smaller than absImaginaryThreshold.
      * absImaginaryThreshold takes the dummy_precision associated
      * with the _Scalar template parameter of the PolynomialSolver class as the default value.
      *
      * \param[out] bi_seq : the back insertion sequence (stl concept)
      * \param[in]  absImaginaryThreshold : the maximum bound of the imaginary part of a complex
      *  number that is considered as real.
      * */
     template<typename Stl_back_insertion_sequence>
     inline void realRoots( Stl_back_insertion_sequence& bi_seq,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       using std::abs;
       bi_seq.clear();
       for(Index i=0; i<m_roots.size(); ++i )
       {
         if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){
           bi_seq.push_back( m_roots[i].real() ); }
       }
     }
 
   protected:
     template<typename squaredNormBinaryPredicate>
     inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
     {
       Index res=0;
       RealScalar norm2 = numext::abs2( m_roots[0] );
       for( Index i=1; i<m_roots.size(); ++i )
       {
         const RealScalar currNorm2 = numext::abs2( m_roots[i] );
         if( pred( currNorm2, norm2 ) ){
           res=i; norm2=currNorm2; }
       }
       return m_roots[res];
     }
 
   public:
     /**
      * \returns the complex root with greatest norm.
      */
     inline const RootType& greatestRoot() const
     {
       std::greater<RealScalar> greater;
       return selectComplexRoot_withRespectToNorm( greater );
     }
 
     /**
      * \returns the complex root with smallest norm.
      */
     inline const RootType& smallestRoot() const
     {
       std::less<RealScalar> less;
       return selectComplexRoot_withRespectToNorm( less );
     }
 
   protected:
     template<typename squaredRealPartBinaryPredicate>
     inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
         squaredRealPartBinaryPredicate& pred,
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       using std::abs;
       hasArealRoot = false;
       Index res=0;
       RealScalar abs2(0);
 
       for( Index i=0; i<m_roots.size(); ++i )
       {
         if( abs( m_roots[i].imag() ) <= absImaginaryThreshold )
         {
           if( !hasArealRoot )
           {
             hasArealRoot = true;
             res = i;
             abs2 = m_roots[i].real() * m_roots[i].real();
           }
           else
           {
             const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
             if( pred( currAbs2, abs2 ) )
             {
               abs2 = currAbs2;
               res = i;
             }
           }
         }
         else if(!hasArealRoot)
         {
           if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
             res = i;}
         }
       }
       return numext::real_ref(m_roots[res]);
     }
 
 
     template<typename RealPartBinaryPredicate>
     inline const RealScalar& selectRealRoot_withRespectToRealPart(
         RealPartBinaryPredicate& pred,
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       using std::abs;
       hasArealRoot = false;
       Index res=0;
       RealScalar val(0);
 
       for( Index i=0; i<m_roots.size(); ++i )
       {
         if( abs( m_roots[i].imag() ) <= absImaginaryThreshold )
         {
           if( !hasArealRoot )
           {
             hasArealRoot = true;
             res = i;
             val = m_roots[i].real();
           }
           else
           {
             const RealScalar curr = m_roots[i].real();
             if( pred( curr, val ) )
             {
               val = curr;
               res = i;
             }
           }
         }
         else
         {
           if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
             res = i; }
         }
       }
       return numext::real_ref(m_roots[res]);
     }
 
   public:
     /**
      * \returns a real root with greatest absolute magnitude.
      * A real root is defined as the real part of a complex root with absolute imaginary
      * part smallest than absImaginaryThreshold.
      * absImaginaryThreshold takes the dummy_precision associated
      * with the _Scalar template parameter of the PolynomialSolver class as the default value.
      * If no real root is found the boolean hasArealRoot is set to false and the real part of
      * the root with smallest absolute imaginary part is returned instead.
      *
      * \param[out] hasArealRoot : boolean true if a real root is found according to the
      *  absImaginaryThreshold criterion, false otherwise.
      * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
      *  whether or not a root is real.
      */
     inline const RealScalar& absGreatestRealRoot(
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       std::greater<RealScalar> greater;
       return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
     }
 
 
     /**
      * \returns a real root with smallest absolute magnitude.
      * A real root is defined as the real part of a complex root with absolute imaginary
      * part smallest than absImaginaryThreshold.
      * absImaginaryThreshold takes the dummy_precision associated
      * with the _Scalar template parameter of the PolynomialSolver class as the default value.
      * If no real root is found the boolean hasArealRoot is set to false and the real part of
      * the root with smallest absolute imaginary part is returned instead.
      *
      * \param[out] hasArealRoot : boolean true if a real root is found according to the
      *  absImaginaryThreshold criterion, false otherwise.
      * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
      *  whether or not a root is real.
      */
     inline const RealScalar& absSmallestRealRoot(
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       std::less<RealScalar> less;
       return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
     }
 
 
     /**
      * \returns the real root with greatest value.
      * A real root is defined as the real part of a complex root with absolute imaginary
      * part smallest than absImaginaryThreshold.
      * absImaginaryThreshold takes the dummy_precision associated
      * with the _Scalar template parameter of the PolynomialSolver class as the default value.
      * If no real root is found the boolean hasArealRoot is set to false and the real part of
      * the root with smallest absolute imaginary part is returned instead.
      *
      * \param[out] hasArealRoot : boolean true if a real root is found according to the
      *  absImaginaryThreshold criterion, false otherwise.
      * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
      *  whether or not a root is real.
      */
     inline const RealScalar& greatestRealRoot(
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       std::greater<RealScalar> greater;
       return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
     }
 
 
     /**
      * \returns the real root with smallest value.
      * A real root is defined as the real part of a complex root with absolute imaginary
      * part smallest than absImaginaryThreshold.
      * absImaginaryThreshold takes the dummy_precision associated
      * with the _Scalar template parameter of the PolynomialSolver class as the default value.
      * If no real root is found the boolean hasArealRoot is set to false and the real part of
      * the root with smallest absolute imaginary part is returned instead.
      *
      * \param[out] hasArealRoot : boolean true if a real root is found according to the
      *  absImaginaryThreshold criterion, false otherwise.
      * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
      *  whether or not a root is real.
      */
     inline const RealScalar& smallestRealRoot(
         bool& hasArealRoot,
         const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
     {
       std::less<RealScalar> less;
       return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
     }
 
   protected:
     RootsType               m_roots;
 };
 
 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE )  \
   typedef typename BASE::Scalar                 Scalar;       \
   typedef typename BASE::RealScalar             RealScalar;   \
   typedef typename BASE::RootType               RootType;     \
   typedef typename BASE::RootsType              RootsType;
 
 
 
 /** \ingroup Polynomials_Module
   *
   * \class PolynomialSolver
   *
   * \brief A polynomial solver
   *
   * Computes the complex roots of a real polynomial.
   *
   * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients
   * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic.
   *             Notice that the number of polynomial coefficients is _Deg+1.
   *
   * This class implements a polynomial solver and provides convenient methods such as
   * - real roots,
   * - greatest, smallest complex roots,
   * - real roots with greatest, smallest absolute real value.
   * - greatest, smallest real roots.
   *
   * WARNING: this polynomial solver is experimental, part of the unsupported Eigen modules.
   *
   *
   * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of
   * the polynomial to compute its roots.
   * This supposes that the complex moduli of the roots are all distinct: e.g. there should
   * be no multiple roots or conjugate roots for instance.
   * With 32bit (float) floating types this problem shows up frequently.
   * However, almost always, correct accuracy is reached even in these cases for 64bit
   * (double) floating types and small polynomial degree (<20).
   */
 template<typename _Scalar, int _Deg>
 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
 {
   public:
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
 
     typedef PolynomialSolverBase<_Scalar,_Deg>    PS_Base;
     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
 
     typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
                                           ComplexEigenSolver<CompanionMatrixType>,
                                           EigenSolver<CompanionMatrixType> >::type EigenSolverType;
     typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar> >::type ComplexScalar;
 
   public:
     /** Computes the complex roots of a new polynomial. */
     template< typename OtherPolynomial >
     void compute( const OtherPolynomial& poly )
     {
       eigen_assert( Scalar(0) != poly[poly.size()-1] );
       eigen_assert( poly.size() > 1 );
       if(poly.size() >  2 )
       {
         internal::companion<Scalar,_Deg> companion( poly );
         companion.balance();
         m_eigenSolver.compute( companion.denseMatrix() );
         m_roots = m_eigenSolver.eigenvalues();
         // cleanup noise in imaginary part of real roots:
         // if the imaginary part is rather small compared to the real part
         // and that cancelling the imaginary part yield a smaller evaluation,
         // then it's safe to keep the real part only.
         RealScalar coarse_prec = RealScalar(std::pow(4,poly.size()+1))*NumTraits<RealScalar>::epsilon();
         for(Index i = 0; i<m_roots.size(); ++i)
         {
           if( internal::isMuchSmallerThan(numext::abs(numext::imag(m_roots[i])),
                                           numext::abs(numext::real(m_roots[i])),
                                           coarse_prec) )
           {
             ComplexScalar as_real_root = ComplexScalar(numext::real(m_roots[i]));
             if(    numext::abs(poly_eval(poly, as_real_root))
                 <= numext::abs(poly_eval(poly, m_roots[i])))
             {
               m_roots[i] = as_real_root;
             }
           }
         }
       }
       else if(poly.size () == 2)
       {
         m_roots.resize(1);
         m_roots[0] = -poly[0]/poly[1];
       }
     }
 
   public:
     template< typename OtherPolynomial >
     inline PolynomialSolver( const OtherPolynomial& poly ){
       compute( poly ); }
 
     inline PolynomialSolver(){}
 
   protected:
     using                   PS_Base::m_roots;
     EigenSolverType         m_eigenSolver;
 };
 
 
 template< typename _Scalar >
 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
 {
   public:
     typedef PolynomialSolverBase<_Scalar,1>    PS_Base;
     EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
 
   public:
     /** Computes the complex roots of a new polynomial. */
     template< typename OtherPolynomial >
     void compute( const OtherPolynomial& poly )
     {
       eigen_assert( poly.size() == 2 );
       eigen_assert( Scalar(0) != poly[1] );
       m_roots[0] = -poly[0]/poly[1];
     }
 
   public:
     template< typename OtherPolynomial >
     inline PolynomialSolver( const OtherPolynomial& poly ){
       compute( poly ); }
 
     inline PolynomialSolver(){}
 
   protected:
     using                   PS_Base::m_roots;
 };
 
 } // end namespace Eigen
 
 #endif // EIGEN_POLYNOMIAL_SOLVER_H

This page was automatically generated by the 2.3.7 LXR engine. The LXR team