PolynomialUtils.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_UTILS_H
26 #define EIGEN_POLYNOMIAL_UTILS_H
27 
28 namespace Eigen {
29 
41 template <typename Polynomials, typename T>
42 inline
43 T poly_eval_horner( const Polynomials& poly, const T& x )
44 {
45  T val=poly[poly.size()-1];
46  for(DenseIndex i=poly.size()-2; i>=0; --i ){
47  val = val*x + poly[i]; }
48  return val;
49 }
50 
59 template <typename Polynomials, typename T>
60 inline
61 T poly_eval( const Polynomials& poly, const T& x )
62 {
63  typedef typename NumTraits<T>::Real Real;
64 
65  if( internal::abs2( x ) <= Real(1) ){
66  return poly_eval_horner( poly, x ); }
67  else
68  {
69  T val=poly[0];
70  T inv_x = T(1)/x;
71  for( DenseIndex i=1; i<poly.size(); ++i ){
72  val = val*inv_x + poly[i]; }
73 
74  return std::pow(x,(T)(poly.size()-1)) * val;
75  }
76 }
77 
88 template <typename Polynomial>
89 inline
90 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly )
91 {
92  typedef typename Polynomial::Scalar Scalar;
93  typedef typename NumTraits<Scalar>::Real Real;
94 
95  assert( Scalar(0) != poly[poly.size()-1] );
96  const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
97  Real cb(0);
98 
99  for( DenseIndex i=0; i<poly.size()-1; ++i ){
100  cb += internal::abs(poly[i]*inv_leading_coeff); }
101  return cb + Real(1);
102 }
103 
110 template <typename Polynomial>
111 inline
112 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly )
113 {
114  typedef typename Polynomial::Scalar Scalar;
115  typedef typename NumTraits<Scalar>::Real Real;
116 
117  DenseIndex i=0;
118  while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
119  if( poly.size()-1 == i ){
120  return Real(1); }
121 
122  const Scalar inv_min_coeff = Scalar(1)/poly[i];
123  Real cb(1);
124  for( DenseIndex j=i+1; j<poly.size(); ++j ){
125  cb += internal::abs(poly[j]*inv_min_coeff); }
126  return Real(1)/cb;
127 }
128 
139 template <typename RootVector, typename Polynomial>
140 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
141 {
142 
143  typedef typename Polynomial::Scalar Scalar;
144 
145  poly.setZero( rv.size()+1 );
146  poly[0] = -rv[0]; poly[1] = Scalar(1);
147  for( DenseIndex i=1; i< rv.size(); ++i )
148  {
149  for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
150  poly[0] = -rv[i]*poly[0];
151  }
152 }
153 
154 } // end namespace Eigen
155 
156 #endif // EIGEN_POLYNOMIAL_UTILS_H