PolynomialUtils.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_UTILS_H
11 #define EIGEN_POLYNOMIAL_UTILS_H
12 
13 namespace Eigen {
14 
26 template <typename Polynomials, typename T>
27 inline
28 T poly_eval_horner( const Polynomials& poly, const T& x )
29 {
30  T val=poly[poly.size()-1];
31  for(DenseIndex i=poly.size()-2; i>=0; --i ){
32  val = val*x + poly[i]; }
33  return val;
34 }
35 
44 template <typename Polynomials, typename T>
45 inline
46 T poly_eval( const Polynomials& poly, const T& x )
47 {
48  typedef typename NumTraits<T>::Real Real;
49 
50  if( numext::abs2( x ) <= Real(1) ){
51  return poly_eval_horner( poly, x ); }
52  else
53  {
54  T val=poly[0];
55  T inv_x = T(1)/x;
56  for( DenseIndex i=1; i<poly.size(); ++i ){
57  val = val*inv_x + poly[i]; }
58 
59  return numext::pow(x,(T)(poly.size()-1)) * val;
60  }
61 }
62 
73 template <typename Polynomial>
74 inline
76 {
77  using std::abs;
78  typedef typename Polynomial::Scalar Scalar;
79  typedef typename NumTraits<Scalar>::Real Real;
80 
81  eigen_assert( Scalar(0) != poly[poly.size()-1] );
82  const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
83  Real cb(0);
84 
85  for( DenseIndex i=0; i<poly.size()-1; ++i ){
86  cb += abs(poly[i]*inv_leading_coeff); }
87  return cb + Real(1);
88 }
89 
96 template <typename Polynomial>
97 inline
99 {
100  using std::abs;
101  typedef typename Polynomial::Scalar Scalar;
102  typedef typename NumTraits<Scalar>::Real Real;
103 
104  DenseIndex i=0;
105  while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
106  if( poly.size()-1 == i ){
107  return Real(1); }
108 
109  const Scalar inv_min_coeff = Scalar(1)/poly[i];
110  Real cb(1);
111  for( DenseIndex j=i+1; j<poly.size(); ++j ){
112  cb += abs(poly[j]*inv_min_coeff); }
113  return Real(1)/cb;
114 }
115 
126 template <typename RootVector, typename Polynomial>
127 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
128 {
129 
130  typedef typename Polynomial::Scalar Scalar;
131 
132  poly.setZero( rv.size()+1 );
133  poly[0] = -rv[0]; poly[1] = Scalar(1);
134  for( DenseIndex i=1; i< rv.size(); ++i )
135  {
136  for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
137  poly[0] = -rv[i]*poly[0];
138  }
139 }
140 
141 } // end namespace Eigen
142 
143 #endif // EIGEN_POLYNOMIAL_UTILS_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
T poly_eval_horner(const Polynomials &poly, const T &x)
T poly_eval(const Polynomials &poly, const T &x)
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
void roots_to_monicPolynomial(const RootVector &rv, Polynomial &poly)
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::Triplet< double > T
EIGEN_DEVICE_FUNC internal::pow_impl< ScalarX, ScalarY >::result_type pow(const ScalarX &x, const ScalarY &y)
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:66
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
NumTraits< typename Polynomial::Scalar >::Real cauchy_min_bound(const Polynomial &poly)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
EIGEN_DEVICE_FUNC bool abs2(bool x)
NumTraits< typename Polynomial::Scalar >::Real cauchy_max_bound(const Polynomial &poly)
std::ptrdiff_t j


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:14