PolynomialSolver.h
Go to the documentation of this file.
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:
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()-1); }
45 
46  public:
47  template< typename OtherPolynomial >
48  inline PolynomialSolverBase( const OtherPolynomial& poly ){
49  setPolynomial( poly() ); }
50 
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  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;
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:
100  inline const RootType& greatestRoot() const
101  {
102  std::greater<Scalar> greater;
103  return selectComplexRoot_withRespectToNorm( greater );
104  }
105 
109  inline const RootType& smallestRoot() const
110  {
111  std::less<Scalar> less;
113  }
114 
115  protected:
116  template<typename squaredRealPartBinaryPredicate>
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>
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:
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 
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 
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 
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:
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 
330 template< typename _Scalar, int _Deg >
331 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
332 {
333  public:
335 
336  typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base;
338 
339  typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
340  typedef EigenSolver<CompanionMatrixType> EigenSolverType;
341 
342  public:
344  template< typename OtherPolynomial >
345  void compute( const OtherPolynomial& poly )
346  {
347  eigen_assert( Scalar(0) != poly[poly.size()-1] );
348  eigen_assert( poly.size() > 1 );
349  if(poly.size() > 2 )
350  {
351  internal::companion<Scalar,_Deg> companion( poly );
352  companion.balance();
353  m_eigenSolver.compute( companion.denseMatrix() );
354  m_roots = m_eigenSolver.eigenvalues();
355  }
356  else if(poly.size () == 2)
357  {
358  m_roots.resize(1);
359  m_roots[0] = -poly[0]/poly[1];
360  }
361  }
362 
363  public:
364  template< typename OtherPolynomial >
365  inline PolynomialSolver( const OtherPolynomial& poly ){
366  compute( poly ); }
367 
368  inline PolynomialSolver(){}
369 
370  protected:
371  using PS_Base::m_roots;
373 };
374 
375 
376 template< typename _Scalar >
377 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
378 {
379  public:
382 
383  public:
385  template< typename OtherPolynomial >
386  void compute( const OtherPolynomial& poly )
387  {
388  eigen_assert( poly.size() == 2 );
389  eigen_assert( Scalar(0) != poly[1] );
390  m_roots[0] = -poly[0]/poly[1];
391  }
392 
393  public:
394  template< typename OtherPolynomial >
395  inline PolynomialSolver( const OtherPolynomial& poly ){
396  compute( poly ); }
397 
398  inline PolynomialSolver(){}
399 
400  protected:
401  using PS_Base::m_roots;
402 };
403 
404 } // end namespace Eigen
405 
406 #endif // EIGEN_POLYNOMIAL_SOLVER_H
EIGEN_DEVICE_FUNC internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
const RealScalar & greatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
PolynomialSolverBase< _Scalar, 1 > PS_Base
DerType::Scalar imag(const AutoDiffScalar< DerType > &)
Definition: Half.h:150
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Defined to be inherited by polynomial solvers: it provides convenient methods such as...
#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar, Size)
Definition: Memory.h:691
PolynomialSolver(const OtherPolynomial &poly)
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
const RootType & smallestRoot() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
const RealScalar & smallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
const RealScalar & absSmallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
double norm2(const Point2 &p, OptionalJacobian< 1, 2 > H)
Distance of the point from the origin, with Jacobian.
Definition: Point2.cpp:27
DenseCompanionMatrixType denseMatrix() const
Definition: Companion.h:88
#define eigen_assert(x)
Definition: Macros.h:579
#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(BASE)
PolynomialSolver(const OtherPolynomial &poly)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
const RootType & greatestRoot() const
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:25
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
void setPolynomial(const OtherPolynomial &poly)
const RealScalar & absGreatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
PolynomialSolverBase(const OtherPolynomial &poly)
const RootsType & roots() const
EigenSolverType m_eigenSolver
const int Dynamic
Definition: Constants.h:21
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
EIGEN_DONT_INLINE void compute(Solver &solver, const MatrixType &A)
The matrix class, also used for vectors and row-vectors.
const RealScalar & selectRealRoot_withRespectToAbsRealPart(squaredRealPartBinaryPredicate &pred, bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
#define abs(x)
Definition: datatypes.h:17
A polynomial solver.
const RealScalar & selectRealRoot_withRespectToRealPart(RealPartBinaryPredicate &pred, bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
const RootType & selectComplexRoot_withRespectToNorm(squaredNormBinaryPredicate &pred) const


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:27