You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
428 lines
15 KiB
428 lines
15 KiB
// 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 the mozilla.org home page |
|
|
|
#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
|
|
|