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()); }
47 template<
typename OtherPolynomial >
49 setPolynomial( poly() ); }
68 template<
typename Stl_back_insertion_sequence>
69 inline void realRoots( Stl_back_insertion_sequence& bi_seq,
73 for(Index i=0; i<m_roots.size(); ++i )
75 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
76 bi_seq.push_back( m_roots[i].real() ); }
81 template<
typename squaredNormBinaryPredicate>
82 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred )
const
85 RealScalar norm2 = internal::abs2( m_roots[0] );
86 for( Index i=1; i<m_roots.size(); ++i )
88 const RealScalar currNorm2 = internal::abs2( m_roots[i] );
89 if( pred( currNorm2, norm2 ) ){
90 res=i; norm2=currNorm2; }
101 std::greater<Scalar> greater;
102 return selectComplexRoot_withRespectToNorm( greater );
110 std::less<Scalar> less;
111 return selectComplexRoot_withRespectToNorm( less );
115 template<
typename squaredRealPartBinaryPredicate>
116 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
117 squaredRealPartBinaryPredicate& pred,
121 hasArealRoot =
false;
125 for( Index i=0; i<m_roots.size(); ++i )
127 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
133 abs2 = m_roots[i].real() * m_roots[i].real();
137 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
138 if( pred( currAbs2, abs2 ) )
147 if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
151 return internal::real_ref(m_roots[res]);
155 template<
typename RealPartBinaryPredicate>
156 inline const RealScalar& selectRealRoot_withRespectToRealPart(
157 RealPartBinaryPredicate& pred,
159 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() )
const
161 hasArealRoot =
false;
165 for( Index i=0; i<m_roots.size(); ++i )
167 if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
173 val = m_roots[i].real();
177 const RealScalar curr = m_roots[i].real();
178 if( pred( curr, val ) )
187 if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
191 return internal::real_ref(m_roots[res]);
213 std::greater<Scalar> greater;
214 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
236 std::less<Scalar> less;
237 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
259 std::greater<Scalar> greater;
260 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
282 std::less<Scalar> less;
283 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
290 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
291 typedef typename BASE::Scalar Scalar; \
292 typedef typename BASE::RealScalar RealScalar; \
293 typedef typename BASE::RootType RootType; \
294 typedef typename BASE::RootsType RootsType;
327 template<
typename _Scalar,
int _Deg >
331 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==
Dynamic ?
Dynamic : _Deg)
334 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(
PS_Base )
341 template<
typename OtherPolynomial >
344 assert( Scalar(0) != poly[poly.size()-1] );
345 internal::companion<Scalar,_Deg> companion( poly );
347 m_eigenSolver.
compute( companion.denseMatrix() );
352 template<
typename OtherPolynomial >
356 inline PolynomialSolver(){}
359 using PS_Base::m_roots;
360 EigenSolverType m_eigenSolver;
364 template<
typename _Scalar >
365 class PolynomialSolver<_Scalar,1> :
public PolynomialSolverBase<_Scalar,1>
368 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
369 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
373 template< typename OtherPolynomial >
374 void compute( const OtherPolynomial& poly )
376 assert( Scalar(0) != poly[poly.size()-1] );
377 m_roots[0] = -poly[0]/poly[poly.size()-1];
381 using PS_Base::m_roots;
386 #endif // EIGEN_POLYNOMIAL_SOLVER_H