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 >
setPolynomial(const OtherPolynomial & poly)43     inline void setPolynomial( const OtherPolynomial& poly ){
44       m_roots.resize(poly.size()); }
45 
46   public:
47     template< typename OtherPolynomial >
PolynomialSolverBase(const OtherPolynomial & poly)48     inline PolynomialSolverBase( const OtherPolynomial& poly ){
49       setPolynomial( poly() ); }
50 
PolynomialSolverBase()51     inline PolynomialSolverBase(){}
52 
53   public:
54     /** \returns the complex roots of the polynomial */
roots()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>
selectComplexRoot_withRespectToNorm(squaredNormBinaryPredicate & pred)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      */
greatestRoot()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      */
smallestRoot()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 >
compute(const OtherPolynomial & poly)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 >
PolynomialSolver(const OtherPolynomial & poly)356     inline PolynomialSolver( const OtherPolynomial& poly ){
357       compute( poly ); }
358 
PolynomialSolver()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;
EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(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
390