$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_POLYNOMIAL_UTILS_H 00011 #define EIGEN_POLYNOMIAL_UTILS_H 00012 00013 namespace Eigen { 00014 00026 template <typename Polynomials, typename T> 00027 inline 00028 T poly_eval_horner( const Polynomials& poly, const T& x ) 00029 { 00030 T val=poly[poly.size()-1]; 00031 for(DenseIndex i=poly.size()-2; i>=0; --i ){ 00032 val = val*x + poly[i]; } 00033 return val; 00034 } 00035 00044 template <typename Polynomials, typename T> 00045 inline 00046 T poly_eval( const Polynomials& poly, const T& x ) 00047 { 00048 typedef typename NumTraits<T>::Real Real; 00049 00050 if( numext::abs2( x ) <= Real(1) ){ 00051 return poly_eval_horner( poly, x ); } 00052 else 00053 { 00054 T val=poly[0]; 00055 T inv_x = T(1)/x; 00056 for( DenseIndex i=1; i<poly.size(); ++i ){ 00057 val = val*inv_x + poly[i]; } 00058 00059 return std::pow(x,(T)(poly.size()-1)) * val; 00060 } 00061 } 00062 00073 template <typename Polynomial> 00074 inline 00075 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly ) 00076 { 00077 using std::abs; 00078 typedef typename Polynomial::Scalar Scalar; 00079 typedef typename NumTraits<Scalar>::Real Real; 00080 00081 eigen_assert( Scalar(0) != poly[poly.size()-1] ); 00082 const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1]; 00083 Real cb(0); 00084 00085 for( DenseIndex i=0; i<poly.size()-1; ++i ){ 00086 cb += abs(poly[i]*inv_leading_coeff); } 00087 return cb + Real(1); 00088 } 00089 00096 template <typename Polynomial> 00097 inline 00098 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly ) 00099 { 00100 using std::abs; 00101 typedef typename Polynomial::Scalar Scalar; 00102 typedef typename NumTraits<Scalar>::Real Real; 00103 00104 DenseIndex i=0; 00105 while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; } 00106 if( poly.size()-1 == i ){ 00107 return Real(1); } 00108 00109 const Scalar inv_min_coeff = Scalar(1)/poly[i]; 00110 Real cb(1); 00111 for( DenseIndex j=i+1; j<poly.size(); ++j ){ 00112 cb += abs(poly[j]*inv_min_coeff); } 00113 return Real(1)/cb; 00114 } 00115 00126 template <typename RootVector, typename Polynomial> 00127 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) 00128 { 00129 00130 typedef typename Polynomial::Scalar Scalar; 00131 00132 poly.setZero( rv.size()+1 ); 00133 poly[0] = -rv[0]; poly[1] = Scalar(1); 00134 for( DenseIndex i=1; i< rv.size(); ++i ) 00135 { 00136 for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; } 00137 poly[0] = -rv[i]*poly[0]; 00138 } 00139 } 00140 00141 } // end namespace Eigen 00142 00143 #endif // EIGEN_POLYNOMIAL_UTILS_H