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 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \ingroup Polynomials_Module
|
---|
16 | * \class PolynomialSolverBase.
|
---|
17 | *
|
---|
18 | * \brief Defined to be inherited by polynomial solvers: it provides
|
---|
19 | * convenient methods such as
|
---|
20 | * - real roots,
|
---|
21 | * - greatest, smallest complex roots,
|
---|
22 | * - real roots with greatest, smallest absolute real value,
|
---|
23 | * - greatest, smallest real roots.
|
---|
24 | *
|
---|
25 | * It stores the set of roots as a vector of complexes.
|
---|
26 | *
|
---|
27 | */
|
---|
28 | template< typename _Scalar, int _Deg >
|
---|
29 | class PolynomialSolverBase
|
---|
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;
|
---|
37 | typedef Matrix<RootType,_Deg,1> RootsType;
|
---|
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()); }
|
---|
45 |
|
---|
46 | public:
|
---|
47 | template< typename OtherPolynomial >
|
---|
48 | inline PolynomialSolverBase( const OtherPolynomial& poly ){
|
---|
49 | setPolynomial( poly() ); }
|
---|
50 |
|
---|
51 | inline PolynomialSolverBase(){}
|
---|
52 |
|
---|
53 | public:
|
---|
54 | /** \returns the complex roots of the polynomial */
|
---|
55 | inline const RootsType& roots() const { return m_roots; }
|
---|
56 |
|
---|
57 | public:
|
---|
58 | /** Clear and fills the back insertion sequence with the real roots of the polynomial
|
---|
59 | * i.e. the real part of the complex roots that have an imaginary part which
|
---|
60 | * absolute value is smaller than absImaginaryThreshold.
|
---|
61 | * absImaginaryThreshold takes the dummy_precision associated
|
---|
62 | * with the _Scalar template parameter of the PolynomialSolver class as the default value.
|
---|
63 | *
|
---|
64 | * \param[out] bi_seq : the back insertion sequence (stl concept)
|
---|
65 | * \param[in] absImaginaryThreshold : the maximum bound of the imaginary part of a complex
|
---|
66 | * number that is considered as real.
|
---|
67 | * */
|
---|
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:
|
---|
97 | /**
|
---|
98 | * \returns the complex root with greatest norm.
|
---|
99 | */
|
---|
100 | inline const RootType& greatestRoot() const
|
---|
101 | {
|
---|
102 | std::greater<Scalar> greater;
|
---|
103 | return selectComplexRoot_withRespectToNorm( greater );
|
---|
104 | }
|
---|
105 |
|
---|
106 | /**
|
---|
107 | * \returns the complex root with smallest norm.
|
---|
108 | */
|
---|
109 | inline const RootType& smallestRoot() const
|
---|
110 | {
|
---|
111 | std::less<Scalar> 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
|
---|
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:
|
---|
198 | /**
|
---|
199 | * \returns a real root with greatest absolute magnitude.
|
---|
200 | * A real root is defined as the real part of a complex root with absolute imaginary
|
---|
201 | * part smallest than absImaginaryThreshold.
|
---|
202 | * absImaginaryThreshold takes the dummy_precision associated
|
---|
203 | * with the _Scalar template parameter of the PolynomialSolver class as the default value.
|
---|
204 | * If no real root is found the boolean hasArealRoot is set to false and the real part of
|
---|
205 | * the root with smallest absolute imaginary part is returned instead.
|
---|
206 | *
|
---|
207 | * \param[out] hasArealRoot : boolean true if a real root is found according to the
|
---|
208 | * absImaginaryThreshold criterion, false otherwise.
|
---|
209 | * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
|
---|
210 | * whether or not a root is real.
|
---|
211 | */
|
---|
212 | inline const RealScalar& absGreatestRealRoot(
|
---|
213 | bool& hasArealRoot,
|
---|
214 | const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
|
---|
215 | {
|
---|
216 | std::greater<Scalar> greater;
|
---|
217 | return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
|
---|
218 | }
|
---|
219 |
|
---|
220 |
|
---|
221 | /**
|
---|
222 | * \returns a real root with smallest absolute magnitude.
|
---|
223 | * A real root is defined as the real part of a complex root with absolute imaginary
|
---|
224 | * part smallest than absImaginaryThreshold.
|
---|
225 | * absImaginaryThreshold takes the dummy_precision associated
|
---|
226 | * with the _Scalar template parameter of the PolynomialSolver class as the default value.
|
---|
227 | * If no real root is found the boolean hasArealRoot is set to false and the real part of
|
---|
228 | * the root with smallest absolute imaginary part is returned instead.
|
---|
229 | *
|
---|
230 | * \param[out] hasArealRoot : boolean true if a real root is found according to the
|
---|
231 | * absImaginaryThreshold criterion, false otherwise.
|
---|
232 | * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
|
---|
233 | * whether or not a root is real.
|
---|
234 | */
|
---|
235 | inline const RealScalar& absSmallestRealRoot(
|
---|
236 | bool& hasArealRoot,
|
---|
237 | const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
|
---|
238 | {
|
---|
239 | std::less<Scalar> less;
|
---|
240 | return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
|
---|
241 | }
|
---|
242 |
|
---|
243 |
|
---|
244 | /**
|
---|
245 | * \returns the real root with greatest value.
|
---|
246 | * A real root is defined as the real part of a complex root with absolute imaginary
|
---|
247 | * part smallest than absImaginaryThreshold.
|
---|
248 | * absImaginaryThreshold takes the dummy_precision associated
|
---|
249 | * with the _Scalar template parameter of the PolynomialSolver class as the default value.
|
---|
250 | * If no real root is found the boolean hasArealRoot is set to false and the real part of
|
---|
251 | * the root with smallest absolute imaginary part is returned instead.
|
---|
252 | *
|
---|
253 | * \param[out] hasArealRoot : boolean true if a real root is found according to the
|
---|
254 | * absImaginaryThreshold criterion, false otherwise.
|
---|
255 | * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
|
---|
256 | * whether or not a root is real.
|
---|
257 | */
|
---|
258 | inline const RealScalar& greatestRealRoot(
|
---|
259 | bool& hasArealRoot,
|
---|
260 | const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
|
---|
261 | {
|
---|
262 | std::greater<Scalar> greater;
|
---|
263 | return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
|
---|
264 | }
|
---|
265 |
|
---|
266 |
|
---|
267 | /**
|
---|
268 | * \returns the real root with smallest value.
|
---|
269 | * A real root is defined as the real part of a complex root with absolute imaginary
|
---|
270 | * part smallest than absImaginaryThreshold.
|
---|
271 | * absImaginaryThreshold takes the dummy_precision associated
|
---|
272 | * with the _Scalar template parameter of the PolynomialSolver class as the default value.
|
---|
273 | * If no real root is found the boolean hasArealRoot is set to false and the real part of
|
---|
274 | * the root with smallest absolute imaginary part is returned instead.
|
---|
275 | *
|
---|
276 | * \param[out] hasArealRoot : boolean true if a real root is found according to the
|
---|
277 | * absImaginaryThreshold criterion, false otherwise.
|
---|
278 | * \param[in] absImaginaryThreshold : threshold on the absolute imaginary part to decide
|
---|
279 | * whether or not a root is real.
|
---|
280 | */
|
---|
281 | inline const RealScalar& smallestRealRoot(
|
---|
282 | bool& hasArealRoot,
|
---|
283 | const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
|
---|
284 | {
|
---|
285 | std::less<Scalar> 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 |
|
---|
301 | /** \ingroup Polynomials_Module
|
---|
302 | *
|
---|
303 | * \class PolynomialSolver
|
---|
304 | *
|
---|
305 | * \brief A polynomial solver
|
---|
306 | *
|
---|
307 | * Computes the complex roots of a real polynomial.
|
---|
308 | *
|
---|
309 | * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients
|
---|
310 | * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic.
|
---|
311 | * Notice that the number of polynomial coefficients is _Deg+1.
|
---|
312 | *
|
---|
313 | * This class implements a polynomial solver and provides convenient methods such as
|
---|
314 | * - real roots,
|
---|
315 | * - greatest, smallest complex roots,
|
---|
316 | * - real roots with greatest, smallest absolute real value.
|
---|
317 | * - greatest, smallest real roots.
|
---|
318 | *
|
---|
319 | * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules.
|
---|
320 | *
|
---|
321 | *
|
---|
322 | * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of
|
---|
323 | * the polynomial to compute its roots.
|
---|
324 | * This supposes that the complex moduli of the roots are all distinct: e.g. there should
|
---|
325 | * be no multiple roots or conjugate roots for instance.
|
---|
326 | * With 32bit (float) floating types this problem shows up frequently.
|
---|
327 | * However, almost always, correct accuracy is reached even in these cases for 64bit
|
---|
328 | * (double) floating types and small polynomial degree (<20).
|
---|
329 | */
|
---|
330 | template< typename _Scalar, int _Deg >
|
---|
331 | class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
|
---|
332 | {
|
---|
333 | public:
|
---|
334 | EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
|
---|
335 |
|
---|
336 | typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base;
|
---|
337 | EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
|
---|
338 |
|
---|
339 | typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
|
---|
340 | typedef EigenSolver<CompanionMatrixType> EigenSolverType;
|
---|
341 |
|
---|
342 | public:
|
---|
343 | /** Computes the complex roots of a new polynomial. */
|
---|
344 | template< typename OtherPolynomial >
|
---|
345 | void compute( const OtherPolynomial& poly )
|
---|
346 | {
|
---|
347 | eigen_assert( Scalar(0) != poly[poly.size()-1] );
|
---|
348 | internal::companion<Scalar,_Deg> companion( poly );
|
---|
349 | companion.balance();
|
---|
350 | m_eigenSolver.compute( companion.denseMatrix() );
|
---|
351 | m_roots = m_eigenSolver.eigenvalues();
|
---|
352 | }
|
---|
353 |
|
---|
354 | public:
|
---|
355 | template< typename OtherPolynomial >
|
---|
356 | inline PolynomialSolver( const OtherPolynomial& poly ){
|
---|
357 | compute( poly ); }
|
---|
358 |
|
---|
359 | inline PolynomialSolver(){}
|
---|
360 |
|
---|
361 | protected:
|
---|
362 | using PS_Base::m_roots;
|
---|
363 | EigenSolverType m_eigenSolver;
|
---|
364 | };
|
---|
365 |
|
---|
366 |
|
---|
367 | template< typename _Scalar >
|
---|
368 | class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
|
---|
369 | {
|
---|
370 | public:
|
---|
371 | typedef PolynomialSolverBase<_Scalar,1> PS_Base;
|
---|
372 | EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
|
---|
373 |
|
---|
374 | public:
|
---|
375 | /** Computes the complex roots of a new polynomial. */
|
---|
376 | template< typename OtherPolynomial >
|
---|
377 | void compute( const OtherPolynomial& poly )
|
---|
378 | {
|
---|
379 | eigen_assert( Scalar(0) != poly[poly.size()-1] );
|
---|
380 | m_roots[0] = -poly[0]/poly[poly.size()-1];
|
---|
381 | }
|
---|
382 |
|
---|
383 | protected:
|
---|
384 | using PS_Base::m_roots;
|
---|
385 | };
|
---|
386 |
|
---|
387 | } // end namespace Eigen
|
---|
388 |
|
---|
389 | #endif // EIGEN_POLYNOMIAL_SOLVER_H
|
---|