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<RealScalar> greater;
103 return selectComplexRoot_withRespectToNorm( greater );
111 std::less<RealScalar> 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 )
135 abs2 = m_roots[i].real() * m_roots[i].real();
139 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
140 if( pred( currAbs2,
abs2 ) )
147 else if(!hasArealRoot)
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<RealScalar> greater;
217 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
239 std::less<RealScalar> less;
240 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
262 std::greater<RealScalar> greater;
263 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
285 std::less<RealScalar> 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 )
340 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
343 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar> >::type ComplexScalar;
347 template<
typename OtherPolynomial >
350 eigen_assert( Scalar(0) != poly[poly.size()-1] );
351 eigen_assert( poly.size() > 1 );
354 internal::companion<Scalar,_Deg> companion( poly );
356 m_eigenSolver.compute( companion.denseMatrix() );
357 m_roots = m_eigenSolver.eigenvalues();
363 for(Index i = 0; i<m_roots.size(); ++i)
365 if( internal::isMuchSmallerThan(numext::abs(numext::imag(m_roots[i])),
366 numext::abs(numext::real(m_roots[i])),
369 ComplexScalar as_real_root = ComplexScalar(numext::real(m_roots[i]));
370 if( numext::abs(
poly_eval(poly, as_real_root))
371 <= numext::abs(
poly_eval(poly, m_roots[i])))
373 m_roots[i] = as_real_root;
378 else if(poly.size () == 2)
381 m_roots[0] = -poly[0]/poly[1];
386 template<
typename OtherPolynomial >
390 inline PolynomialSolver(){}
393 using PS_Base::m_roots;
394 EigenSolverType m_eigenSolver;
398 template<
typename _Scalar >
399 class PolynomialSolver<_Scalar,1> :
public PolynomialSolverBase<_Scalar,1>
402 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
403 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
407 template<
typename OtherPolynomial >
408 void compute(
const OtherPolynomial& poly )
410 eigen_assert( poly.size() == 2 );
411 eigen_assert( Scalar(0) != poly[1] );
412 m_roots[0] = -poly[0]/poly[1];
416 template<
typename OtherPolynomial >
417 inline PolynomialSolver(
const OtherPolynomial& poly ){
420 inline PolynomialSolver(){}
423 using PS_Base::m_roots;
428 #endif // EIGEN_POLYNOMIAL_SOLVER_H