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
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 */
28template< typename _Scalar, int _Deg >
29class 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  */
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
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
367template< typename _Scalar >
368class 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
390