Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
PolynomialSolver.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_POLYNOMIAL_SOLVER_H
11#define EIGEN_POLYNOMIAL_SOLVER_H
12
13namespace Eigen {
14
28template< typename _Scalar, int _Deg >
30{
31 public:
32 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
33
34 typedef _Scalar Scalar;
35 typedef typename NumTraits<Scalar>::Real RealScalar;
36 typedef std::complex<RealScalar> RootType;
38
39 typedef DenseIndex Index;
40
41 protected:
42 template< typename OtherPolynomial >
43 inline void setPolynomial( const OtherPolynomial& poly ){
44 m_roots.resize(poly.size()-1); }
45
46 public:
47 template< typename OtherPolynomial >
48 inline PolynomialSolverBase( const OtherPolynomial& poly ){
49 setPolynomial( poly() ); }
50
51 inline PolynomialSolverBase(){}
52
53 public:
55 inline const RootsType& roots() const { return m_roots; }
56
57 public:
68 template<typename Stl_back_insertion_sequence>
69 inline void realRoots( Stl_back_insertion_sequence& bi_seq,
70 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
71 {
72 using std::abs;
73 bi_seq.clear();
74 for(Index i=0; i<m_roots.size(); ++i )
75 {
76 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){
77 bi_seq.push_back( m_roots[i].real() ); }
78 }
79 }
80
81 protected:
82 template<typename squaredNormBinaryPredicate>
83 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
84 {
85 Index res=0;
86 RealScalar norm2 = numext::abs2( m_roots[0] );
87 for( Index i=1; i<m_roots.size(); ++i )
88 {
89 const RealScalar currNorm2 = numext::abs2( m_roots[i] );
90 if( pred( currNorm2, norm2 ) ){
91 res=i; norm2=currNorm2; }
92 }
93 return m_roots[res];
94 }
95
96 public:
100 inline const RootType& greatestRoot() const
101 {
102 std::greater<RealScalar> greater;
103 return selectComplexRoot_withRespectToNorm( greater );
104 }
105
109 inline const RootType& smallestRoot() const
110 {
111 std::less<RealScalar> less;
112 return selectComplexRoot_withRespectToNorm( less );
113 }
114
115 protected:
116 template<typename squaredRealPartBinaryPredicate>
117 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
118 squaredRealPartBinaryPredicate& pred,
119 bool& hasArealRoot,
120 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
121 {
122 using std::abs;
123 hasArealRoot = false;
124 Index res=0;
125 RealScalar abs2(0);
126
127 for( Index i=0; i<m_roots.size(); ++i )
128 {
129 if( abs( m_roots[i].imag() ) <= absImaginaryThreshold )
130 {
131 if( !hasArealRoot )
132 {
133 hasArealRoot = true;
134 res = i;
135 abs2 = m_roots[i].real() * m_roots[i].real();
136 }
137 else
138 {
139 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
140 if( pred( currAbs2, abs2 ) )
141 {
142 abs2 = currAbs2;
143 res = i;
144 }
145 }
146 }
147 else if(!hasArealRoot)
148 {
149 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
150 res = i;}
151 }
152 }
153 return numext::real_ref(m_roots[res]);
154 }
155
156
157 template<typename RealPartBinaryPredicate>
158 inline const RealScalar& selectRealRoot_withRespectToRealPart(
159 RealPartBinaryPredicate& pred,
160 bool& hasArealRoot,
161 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
162 {
163 using std::abs;
164 hasArealRoot = false;
165 Index res=0;
166 RealScalar val(0);
167
168 for( Index i=0; i<m_roots.size(); ++i )
169 {
170 if( abs( m_roots[i].imag() ) <= absImaginaryThreshold )
171 {
172 if( !hasArealRoot )
173 {
174 hasArealRoot = true;
175 res = i;
176 val = m_roots[i].real();
177 }
178 else
179 {
180 const RealScalar curr = m_roots[i].real();
181 if( pred( curr, val ) )
182 {
183 val = curr;
184 res = i;
185 }
186 }
187 }
188 else
189 {
190 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
191 res = i; }
192 }
193 }
194 return numext::real_ref(m_roots[res]);
195 }
196
197 public:
212 inline const RealScalar& absGreatestRealRoot(
213 bool& hasArealRoot,
214 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
215 {
216 std::greater<RealScalar> greater;
217 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
218 }
219
220
235 inline const RealScalar& absSmallestRealRoot(
236 bool& hasArealRoot,
237 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
238 {
239 std::less<RealScalar> less;
240 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
241 }
242
243
258 inline const RealScalar& greatestRealRoot(
259 bool& hasArealRoot,
260 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
261 {
262 std::greater<RealScalar> greater;
263 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
264 }
265
266
281 inline const RealScalar& smallestRealRoot(
282 bool& hasArealRoot,
283 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
284 {
285 std::less<RealScalar> less;
286 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
287 }
288
289 protected:
290 RootsType m_roots;
291};
292
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;
298
299
300
330template<typename _Scalar, int _Deg>
331class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
332{
333 public:
334 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
335
337 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
338
340 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
342 EigenSolver<CompanionMatrixType> >::type EigenSolverType;
343 typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar> >::type ComplexScalar;
344
345 public:
347 template< typename OtherPolynomial >
348 void compute( const OtherPolynomial& poly )
349 {
350 eigen_assert( Scalar(0) != poly[poly.size()-1] );
351 eigen_assert( poly.size() > 1 );
352 if(poly.size() > 2 )
353 {
354 internal::companion<Scalar,_Deg> companion( poly );
355 companion.balance();
356 m_eigenSolver.compute( companion.denseMatrix() );
357 m_roots = m_eigenSolver.eigenvalues();
358 // cleanup noise in imaginary part of real roots:
359 // if the imaginary part is rather small compared to the real part
360 // and that cancelling the imaginary part yield a smaller evaluation,
361 // then it's safe to keep the real part only.
362 RealScalar coarse_prec = RealScalar(std::pow(4,poly.size()+1))*NumTraits<RealScalar>::epsilon();
363 for(Index i = 0; i<m_roots.size(); ++i)
364 {
365 if( internal::isMuchSmallerThan(numext::abs(numext::imag(m_roots[i])),
366 numext::abs(numext::real(m_roots[i])),
367 coarse_prec) )
368 {
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])))
372 {
373 m_roots[i] = as_real_root;
374 }
375 }
376 }
377 }
378 else if(poly.size () == 2)
379 {
380 m_roots.resize(1);
381 m_roots[0] = -poly[0]/poly[1];
382 }
383 }
384
385 public:
386 template< typename OtherPolynomial >
387 inline PolynomialSolver( const OtherPolynomial& poly ){
388 compute( poly ); }
389
390 inline PolynomialSolver(){}
391
392 protected:
393 using PS_Base::m_roots;
394 EigenSolverType m_eigenSolver;
395};
396
397
398template< typename _Scalar >
399class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
400{
401 public:
402 typedef PolynomialSolverBase<_Scalar,1> PS_Base;
403 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
404
405 public:
407 template< typename OtherPolynomial >
408 void compute( const OtherPolynomial& poly )
409 {
410 eigen_assert( poly.size() == 2 );
411 eigen_assert( Scalar(0) != poly[1] );
412 m_roots[0] = -poly[0]/poly[1];
413 }
414
415 public:
416 template< typename OtherPolynomial >
417 inline PolynomialSolver( const OtherPolynomial& poly ){
418 compute( poly ); }
419
420 inline PolynomialSolver(){}
421
422 protected:
423 using PS_Base::m_roots;
424};
425
426} // end namespace Eigen
427
428#endif // EIGEN_POLYNOMIAL_SOLVER_H
void resize(Index rows, Index cols)
Defined to be inherited by polynomial solvers: it provides convenient methods such as.
Definition: PolynomialSolver.h:30
const RealScalar & greatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:258
const RealScalar & absSmallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:235
const RootType & smallestRoot() const
Definition: PolynomialSolver.h:109
const RootsType & roots() const
Definition: PolynomialSolver.h:55
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:69
const RealScalar & absGreatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:212
const RootType & greatestRoot() const
Definition: PolynomialSolver.h:100
const RealScalar & smallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Definition: PolynomialSolver.h:281
A polynomial solver.
Definition: PolynomialSolver.h:332
void compute(const OtherPolynomial &poly)
Definition: PolynomialSolver.h:348
T poly_eval(const Polynomials &poly, const T &x)
Definition: PolynomialUtils.h:46
Namespace containing all symbols from the Eigen library.
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs2_op< typename Derived::Scalar >, const Derived > abs2(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_imag_op< typename Derived::Scalar >, const Derived > imag(const Eigen::ArrayBase< Derived > &x)
const int Dynamic