25 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
26 #define EIGEN_POLYNOMIAL_SOLVER_H
43 template<
typename _Scalar,
int _Deg >
47 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==
Dynamic ?
Dynamic : _Deg)
49 typedef _Scalar Scalar;
50 typedef typename NumTraits<Scalar>::Real RealScalar;
51 typedef std::complex<RealScalar> RootType;
52 typedef Matrix<RootType,_Deg,1> RootsType;
54 typedef DenseIndex Index;
57 template<
typename OtherPolynomial >
58 inline void setPolynomial(
const OtherPolynomial& poly ){
59 m_roots.resize(poly.size()); }
62 template<
typename OtherPolynomial >
64 setPolynomial( poly() ); }
70 inline const RootsType&
roots()
const {
return m_roots; }
83 template<
typename Stl_back_insertion_sequence>
84 inline void realRoots( Stl_back_insertion_sequence& bi_seq,
85 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
88 for(Index i=0; i<m_roots.size(); ++i )
90 if( internal::abs( m_roots[i].
imag() ) < absImaginaryThreshold ){
91 bi_seq.push_back( m_roots[i].
real() ); }
96 template<
typename squaredNormBinaryPredicate>
97 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred )
const
100 RealScalar norm2 = internal::abs2( m_roots[0] );
101 for( Index i=1; i<m_roots.size(); ++i )
103 const RealScalar currNorm2 = internal::abs2( m_roots[i] );
104 if( pred( currNorm2, norm2 ) ){
105 res=i; norm2=currNorm2; }
116 std::greater<Scalar> greater;
117 return selectComplexRoot_withRespectToNorm( greater );
125 std::less<Scalar> less;
126 return selectComplexRoot_withRespectToNorm( less );
130 template<
typename squaredRealPartBinaryPredicate>
131 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
132 squaredRealPartBinaryPredicate& pred,
134 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
136 hasArealRoot =
false;
140 for( Index i=0; i<m_roots.size(); ++i )
142 if( internal::abs( m_roots[i].
imag() ) < absImaginaryThreshold )
148 abs2 = m_roots[i].real() * m_roots[i].real();
152 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
153 if( pred( currAbs2, abs2 ) )
162 if( internal::abs( m_roots[i].
imag() ) < internal::abs( m_roots[res].
imag() ) ){
166 return internal::real_ref(m_roots[res]);
170 template<
typename RealPartBinaryPredicate>
171 inline const RealScalar& selectRealRoot_withRespectToRealPart(
172 RealPartBinaryPredicate& pred,
174 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
176 hasArealRoot =
false;
180 for( Index i=0; i<m_roots.size(); ++i )
182 if( internal::abs( m_roots[i].
imag() ) < absImaginaryThreshold )
188 val = m_roots[i].real();
192 const RealScalar curr = m_roots[i].real();
193 if( pred( curr, val ) )
202 if( internal::abs( m_roots[i].
imag() ) < internal::abs( m_roots[res].
imag() ) ){
206 return internal::real_ref(m_roots[res]);
226 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
228 std::greater<Scalar> greater;
229 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
249 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
251 std::less<Scalar> less;
252 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
272 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
274 std::greater<Scalar> greater;
275 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
295 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
297 std::less<Scalar> less;
298 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
305 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
306 typedef typename BASE::Scalar Scalar; \
307 typedef typename BASE::RealScalar RealScalar; \
308 typedef typename BASE::RootType RootType; \
309 typedef typename BASE::RootsType RootsType;
342 template<
typename _Scalar,
int _Deg >
346 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==
Dynamic ?
Dynamic : _Deg)
349 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(
PS_Base )
351 typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
352 typedef EigenSolver<CompanionMatrixType> EigenSolverType;
356 template<
typename OtherPolynomial >
359 assert( Scalar(0) != poly[poly.size()-1] );
360 internal::companion<Scalar,_Deg> companion( poly );
362 m_eigenSolver.compute( companion.denseMatrix() );
363 m_roots = m_eigenSolver.eigenvalues();
367 template<
typename OtherPolynomial >
371 inline PolynomialSolver(){}
374 using PS_Base::m_roots;
375 EigenSolverType m_eigenSolver;
379 template<
typename _Scalar >
380 class PolynomialSolver<_Scalar,1> :
public PolynomialSolverBase<_Scalar,1>
383 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
384 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
388 template< typename OtherPolynomial >
389 void compute( const OtherPolynomial& poly )
391 assert( Scalar(0) != poly[poly.size()-1] );
392 m_roots[0] = -poly[0]/poly[poly.size()-1];
396 using PS_Base::m_roots;
401 #endif // EIGEN_POLYNOMIAL_SOLVER_H