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.
406 lines
14 KiB
406 lines
14 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 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<Scalar> greater;
|
|
return selectComplexRoot_withRespectToNorm( greater );
|
|
}
|
|
|
|
/**
|
|
* \returns the complex root with smallest norm.
|
|
*/
|
|
inline const RootType& smallestRoot() const
|
|
{
|
|
std::less<Scalar> 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( 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<Scalar> 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<Scalar> 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<Scalar> 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<Scalar> 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 EigenSolver<CompanionMatrixType> EigenSolverType;
|
|
|
|
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();
|
|
}
|
|
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
|
|
|