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 
13 namespace Eigen {
14 
28 template< 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()); }
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  bi_seq.clear();
73  for(Index i=0; i<m_roots.size(); ++i )
74  {
75  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
76  bi_seq.push_back( m_roots[i].real() ); }
77  }
78  }
79 
80  protected:
81  template<typename squaredNormBinaryPredicate>
82  inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
83  {
84  Index res=0;
85  RealScalar norm2 = internal::abs2( m_roots[0] );
86  for( Index i=1; i<m_roots.size(); ++i )
87  {
88  const RealScalar currNorm2 = internal::abs2( m_roots[i] );
89  if( pred( currNorm2, norm2 ) ){
90  res=i; norm2=currNorm2; }
91  }
92  return m_roots[res];
93  }
94 
95  public:
99  inline const RootType& greatestRoot() const
100  {
101  std::greater<Scalar> greater;
102  return selectComplexRoot_withRespectToNorm( greater );
103  }
104 
108  inline const RootType& smallestRoot() const
109  {
110  std::less<Scalar> less;
111  return selectComplexRoot_withRespectToNorm( less );
112  }
113 
114  protected:
115  template<typename squaredRealPartBinaryPredicate>
116  inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
117  squaredRealPartBinaryPredicate& pred,
118  bool& hasArealRoot,
119  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
120  {
121  hasArealRoot = false;
122  Index res=0;
123  RealScalar abs2(0);
124 
125  for( Index i=0; i<m_roots.size(); ++i )
126  {
127  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
128  {
129  if( !hasArealRoot )
130  {
131  hasArealRoot = true;
132  res = i;
133  abs2 = m_roots[i].real() * m_roots[i].real();
134  }
135  else
136  {
137  const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
138  if( pred( currAbs2, abs2 ) )
139  {
140  abs2 = currAbs2;
141  res = i;
142  }
143  }
144  }
145  else
146  {
147  if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
148  res = i; }
149  }
150  }
151  return internal::real_ref(m_roots[res]);
152  }
153 
154 
155  template<typename RealPartBinaryPredicate>
156  inline const RealScalar& selectRealRoot_withRespectToRealPart(
157  RealPartBinaryPredicate& pred,
158  bool& hasArealRoot,
159  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
160  {
161  hasArealRoot = false;
162  Index res=0;
163  RealScalar val(0);
164 
165  for( Index i=0; i<m_roots.size(); ++i )
166  {
167  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
168  {
169  if( !hasArealRoot )
170  {
171  hasArealRoot = true;
172  res = i;
173  val = m_roots[i].real();
174  }
175  else
176  {
177  const RealScalar curr = m_roots[i].real();
178  if( pred( curr, val ) )
179  {
180  val = curr;
181  res = i;
182  }
183  }
184  }
185  else
186  {
187  if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
188  res = i; }
189  }
190  }
191  return internal::real_ref(m_roots[res]);
192  }
193 
194  public:
209  inline const RealScalar& absGreatestRealRoot(
210  bool& hasArealRoot,
211  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
212  {
213  std::greater<Scalar> greater;
214  return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
215  }
216 
217 
232  inline const RealScalar& absSmallestRealRoot(
233  bool& hasArealRoot,
234  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
235  {
236  std::less<Scalar> less;
237  return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
238  }
239 
240 
255  inline const RealScalar& greatestRealRoot(
256  bool& hasArealRoot,
257  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
258  {
259  std::greater<Scalar> greater;
260  return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
261  }
262 
263 
278  inline const RealScalar& smallestRealRoot(
279  bool& hasArealRoot,
280  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
281  {
282  std::less<Scalar> less;
283  return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
284  }
285 
286  protected:
287  RootsType m_roots;
288 };
289 
290 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
291  typedef typename BASE::Scalar Scalar; \
292  typedef typename BASE::RealScalar RealScalar; \
293  typedef typename BASE::RootType RootType; \
294  typedef typename BASE::RootsType RootsType;
295 
296 
297 
327 template< typename _Scalar, int _Deg >
328 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
329 {
330  public:
331  EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
332 
334  EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
335 
338 
339  public:
341  template< typename OtherPolynomial >
342  void compute( const OtherPolynomial& poly )
343  {
344  assert( Scalar(0) != poly[poly.size()-1] );
345  internal::companion<Scalar,_Deg> companion( poly );
346  companion.balance();
347  m_eigenSolver.compute( companion.denseMatrix() );
348  m_roots = m_eigenSolver.eigenvalues();
349  }
350 
351  public:
352  template< typename OtherPolynomial >
353  inline PolynomialSolver( const OtherPolynomial& poly ){
354  compute( poly ); }
355 
356  inline PolynomialSolver(){}
357 
358  protected:
359  using PS_Base::m_roots;
360  EigenSolverType m_eigenSolver;
361 };
362 
363 
364 template< typename _Scalar >
365 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
366 {
367  public:
368  typedef PolynomialSolverBase<_Scalar,1> PS_Base;
369  EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
370 
371  public:
373  template< typename OtherPolynomial >
374  void compute( const OtherPolynomial& poly )
375  {
376  assert( Scalar(0) != poly[poly.size()-1] );
377  m_roots[0] = -poly[0]/poly[poly.size()-1];
378  }
379 
380  protected:
381  using PS_Base::m_roots;
382 };
383 
384 } // end namespace Eigen
385 
386 #endif // EIGEN_POLYNOMIAL_SOLVER_H