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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
26 #define EIGEN_POLYNOMIAL_SOLVER_H
27 
28 namespace Eigen {
29 
43 template< typename _Scalar, int _Deg >
45 {
46  public:
47  EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
48 
49  typedef _Scalar Scalar;
50  typedef typename NumTraits<Scalar>::Real RealScalar;
51  typedef std::complex<RealScalar> RootType;
52  typedef Matrix<RootType,_Deg,1> RootsType;
53 
54  typedef DenseIndex Index;
55 
56  protected:
57  template< typename OtherPolynomial >
58  inline void setPolynomial( const OtherPolynomial& poly ){
59  m_roots.resize(poly.size()); }
60 
61  public:
62  template< typename OtherPolynomial >
63  inline PolynomialSolverBase( const OtherPolynomial& poly ){
64  setPolynomial( poly() ); }
65 
66  inline PolynomialSolverBase(){}
67 
68  public:
70  inline const RootsType& roots() const { return m_roots; }
71 
72  public:
83  template<typename Stl_back_insertion_sequence>
84  inline void realRoots( Stl_back_insertion_sequence& bi_seq,
85  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
86  {
87  bi_seq.clear();
88  for(Index i=0; i<m_roots.size(); ++i )
89  {
90  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold ){
91  bi_seq.push_back( m_roots[i].real() ); }
92  }
93  }
94 
95  protected:
96  template<typename squaredNormBinaryPredicate>
97  inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
98  {
99  Index res=0;
100  RealScalar norm2 = internal::abs2( m_roots[0] );
101  for( Index i=1; i<m_roots.size(); ++i )
102  {
103  const RealScalar currNorm2 = internal::abs2( m_roots[i] );
104  if( pred( currNorm2, norm2 ) ){
105  res=i; norm2=currNorm2; }
106  }
107  return m_roots[res];
108  }
109 
110  public:
114  inline const RootType& greatestRoot() const
115  {
116  std::greater<Scalar> greater;
117  return selectComplexRoot_withRespectToNorm( greater );
118  }
119 
123  inline const RootType& smallestRoot() const
124  {
125  std::less<Scalar> less;
126  return selectComplexRoot_withRespectToNorm( less );
127  }
128 
129  protected:
130  template<typename squaredRealPartBinaryPredicate>
131  inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
132  squaredRealPartBinaryPredicate& pred,
133  bool& hasArealRoot,
134  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
135  {
136  hasArealRoot = false;
137  Index res=0;
138  RealScalar abs2(0);
139 
140  for( Index i=0; i<m_roots.size(); ++i )
141  {
142  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
143  {
144  if( !hasArealRoot )
145  {
146  hasArealRoot = true;
147  res = i;
148  abs2 = m_roots[i].real() * m_roots[i].real();
149  }
150  else
151  {
152  const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
153  if( pred( currAbs2, abs2 ) )
154  {
155  abs2 = currAbs2;
156  res = i;
157  }
158  }
159  }
160  else
161  {
162  if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
163  res = i; }
164  }
165  }
166  return internal::real_ref(m_roots[res]);
167  }
168 
169 
170  template<typename RealPartBinaryPredicate>
171  inline const RealScalar& selectRealRoot_withRespectToRealPart(
172  RealPartBinaryPredicate& pred,
173  bool& hasArealRoot,
174  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
175  {
176  hasArealRoot = false;
177  Index res=0;
178  RealScalar val(0);
179 
180  for( Index i=0; i<m_roots.size(); ++i )
181  {
182  if( internal::abs( m_roots[i].imag() ) < absImaginaryThreshold )
183  {
184  if( !hasArealRoot )
185  {
186  hasArealRoot = true;
187  res = i;
188  val = m_roots[i].real();
189  }
190  else
191  {
192  const RealScalar curr = m_roots[i].real();
193  if( pred( curr, val ) )
194  {
195  val = curr;
196  res = i;
197  }
198  }
199  }
200  else
201  {
202  if( internal::abs( m_roots[i].imag() ) < internal::abs( m_roots[res].imag() ) ){
203  res = i; }
204  }
205  }
206  return internal::real_ref(m_roots[res]);
207  }
208 
209  public:
224  inline const RealScalar& absGreatestRealRoot(
225  bool& hasArealRoot,
226  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
227  {
228  std::greater<Scalar> greater;
229  return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
230  }
231 
232 
247  inline const RealScalar& absSmallestRealRoot(
248  bool& hasArealRoot,
249  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
250  {
251  std::less<Scalar> less;
252  return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
253  }
254 
255 
270  inline const RealScalar& greatestRealRoot(
271  bool& hasArealRoot,
272  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
273  {
274  std::greater<Scalar> greater;
275  return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
276  }
277 
278 
293  inline const RealScalar& smallestRealRoot(
294  bool& hasArealRoot,
295  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
296  {
297  std::less<Scalar> less;
298  return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
299  }
300 
301  protected:
302  RootsType m_roots;
303 };
304 
305 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
306  typedef typename BASE::Scalar Scalar; \
307  typedef typename BASE::RealScalar RealScalar; \
308  typedef typename BASE::RootType RootType; \
309  typedef typename BASE::RootsType RootsType;
310 
311 
312 
342 template< typename _Scalar, int _Deg >
343 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
344 {
345  public:
346  EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
347 
349  EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
350 
351  typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
352  typedef EigenSolver<CompanionMatrixType> EigenSolverType;
353 
354  public:
356  template< typename OtherPolynomial >
357  void compute( const OtherPolynomial& poly )
358  {
359  assert( Scalar(0) != poly[poly.size()-1] );
360  internal::companion<Scalar,_Deg> companion( poly );
361  companion.balance();
362  m_eigenSolver.compute( companion.denseMatrix() );
363  m_roots = m_eigenSolver.eigenvalues();
364  }
365 
366  public:
367  template< typename OtherPolynomial >
368  inline PolynomialSolver( const OtherPolynomial& poly ){
369  compute( poly ); }
370 
371  inline PolynomialSolver(){}
372 
373  protected:
374  using PS_Base::m_roots;
375  EigenSolverType m_eigenSolver;
376 };
377 
378 
379 template< typename _Scalar >
380 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
381 {
382  public:
383  typedef PolynomialSolverBase<_Scalar,1> PS_Base;
384  EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
385 
386  public:
388  template< typename OtherPolynomial >
389  void compute( const OtherPolynomial& poly )
390  {
391  assert( Scalar(0) != poly[poly.size()-1] );
392  m_roots[0] = -poly[0]/poly[poly.size()-1];
393  }
394 
395  protected:
396  using PS_Base::m_roots;
397 };
398 
399 } // end namespace Eigen
400 
401 #endif // EIGEN_POLYNOMIAL_SOLVER_H