10 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
11 #define EIGEN_POLYNOMIAL_SOLVER_H
28 template<
typename _Scalar,
int _Deg >
32 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==
Dynamic ?
Dynamic : _Deg)
34 typedef _Scalar Scalar;
36 typedef std::complex<RealScalar> RootType;
39 typedef DenseIndex Index;
42 template<
typename OtherPolynomial >
43 inline void setPolynomial(
const OtherPolynomial& poly ){
44 m_roots.
resize(poly.size()-1); }
47 template<
typename OtherPolynomial >
49 setPolynomial( poly() ); }
68 template<
typename Stl_back_insertion_sequence>
69 inline void realRoots( Stl_back_insertion_sequence& bi_seq,
74 for(Index i=0; i<m_roots.size(); ++i )
76 if(
abs( m_roots[i].
imag() ) < absImaginaryThreshold ){
77 bi_seq.push_back( m_roots[i].
real() ); }
82 template<
typename squaredNormBinaryPredicate>
83 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred )
const
86 RealScalar norm2 = numext::abs2( m_roots[0] );
87 for(
Index i=1; i<m_roots.size(); ++i )
89 const RealScalar currNorm2 = numext::abs2( m_roots[i] );
90 if( pred( currNorm2, norm2 ) ){
91 res=i; norm2=currNorm2; }
102 std::greater<Scalar> greater;
103 return selectComplexRoot_withRespectToNorm( greater );
111 std::less<Scalar> less;
112 return selectComplexRoot_withRespectToNorm( less );
116 template<
typename squaredRealPartBinaryPredicate>
117 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
118 squaredRealPartBinaryPredicate& pred,
123 hasArealRoot =
false;
127 for(
Index i=0; i<m_roots.size(); ++i )
129 if(
abs( m_roots[i].
imag() ) < absImaginaryThreshold )
139 const RealScalar currAbs2 = m_roots[i].
real() * m_roots[i].
real();
140 if( pred( currAbs2,
abs2 ) )
153 return numext::real_ref(m_roots[res]);
157 template<
typename RealPartBinaryPredicate>
158 inline const RealScalar& selectRealRoot_withRespectToRealPart(
159 RealPartBinaryPredicate& pred,
161 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
164 hasArealRoot =
false;
168 for( Index i=0; i<m_roots.size(); ++i )
170 if(
abs( m_roots[i].
imag() ) < absImaginaryThreshold )
176 val = m_roots[i].
real();
180 const RealScalar curr = m_roots[i].
real();
181 if( pred( curr, val ) )
194 return numext::real_ref(m_roots[res]);
216 std::greater<Scalar> greater;
217 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
239 std::less<Scalar> less;
240 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
262 std::greater<Scalar> greater;
263 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
285 std::less<Scalar> less;
286 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
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;
330 template<
typename _Scalar,
int _Deg >
334 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==
Dynamic ?
Dynamic : _Deg)
337 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(
PS_Base )
344 template<
typename OtherPolynomial >
347 eigen_assert( Scalar(0) != poly[poly.size()-1] );
348 eigen_assert( poly.size() > 1 );
351 internal::companion<Scalar,_Deg> companion( poly );
353 m_eigenSolver.
compute( companion.denseMatrix() );
356 else if(poly.size () == 2)
359 m_roots[0] = -poly[0]/poly[1];
364 template<
typename OtherPolynomial >
368 inline PolynomialSolver(){}
371 using PS_Base::m_roots;
372 EigenSolverType m_eigenSolver;
376 template<
typename _Scalar >
377 class PolynomialSolver<_Scalar,1> :
public PolynomialSolverBase<_Scalar,1>
380 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
381 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
385 template< typename OtherPolynomial >
386 void compute( const OtherPolynomial& poly )
388 eigen_assert( poly.size() == 2 );
389 eigen_assert( Scalar(0) != poly[1] );
390 m_roots[0] = -poly[0]/poly[1];
394 template<
typename OtherPolynomial >
395 inline PolynomialSolver(
const OtherPolynomial& poly ){
398 inline PolynomialSolver(){}
401 using PS_Base::m_roots;
406 #endif // EIGEN_POLYNOMIAL_SOLVER_H
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
const RealScalar & greatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:258
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:69
Defined to be inherited by polynomial solvers: it provides convenient methods such as...
Definition: PolynomialSolver.h:29
void resize(Index rows, Index cols)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs2_op< typename Derived::Scalar >, const Derived > abs2(const Eigen::ArrayBase< Derived > &x)
const RootType & smallestRoot() const
Definition: PolynomialSolver.h:109
const RealScalar & smallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:281
const RealScalar & absSmallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:235
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_imag_op< typename Derived::Scalar >, const Derived > imag(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
RealReturnType real() const
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const RootType & greatestRoot() const
Definition: PolynomialSolver.h:100
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const RealScalar & absGreatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:212
const RootsType & roots() const
Definition: PolynomialSolver.h:55
void compute(const OtherPolynomial &poly)
Definition: PolynomialSolver.h:345
A polynomial solver.
Definition: PolynomialSolver.h:331
const EigenvalueType & eigenvalues() const