00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_POLYNOMIAL_UTILS_H
00026 #define EIGEN_POLYNOMIAL_UTILS_H
00027
00028 namespace Eigen {
00029
00041 template <typename Polynomials, typename T>
00042 inline
00043 T poly_eval_horner( const Polynomials& poly, const T& x )
00044 {
00045 T val=poly[poly.size()-1];
00046 for(DenseIndex i=poly.size()-2; i>=0; --i ){
00047 val = val*x + poly[i]; }
00048 return val;
00049 }
00050
00059 template <typename Polynomials, typename T>
00060 inline
00061 T poly_eval( const Polynomials& poly, const T& x )
00062 {
00063 typedef typename NumTraits<T>::Real Real;
00064
00065 if( internal::abs2( x ) <= Real(1) ){
00066 return poly_eval_horner( poly, x ); }
00067 else
00068 {
00069 T val=poly[0];
00070 T inv_x = T(1)/x;
00071 for( DenseIndex i=1; i<poly.size(); ++i ){
00072 val = val*inv_x + poly[i]; }
00073
00074 return std::pow(x,(T)(poly.size()-1)) * val;
00075 }
00076 }
00077
00088 template <typename Polynomial>
00089 inline
00090 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly )
00091 {
00092 typedef typename Polynomial::Scalar Scalar;
00093 typedef typename NumTraits<Scalar>::Real Real;
00094
00095 assert( Scalar(0) != poly[poly.size()-1] );
00096 const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
00097 Real cb(0);
00098
00099 for( DenseIndex i=0; i<poly.size()-1; ++i ){
00100 cb += internal::abs(poly[i]*inv_leading_coeff); }
00101 return cb + Real(1);
00102 }
00103
00110 template <typename Polynomial>
00111 inline
00112 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly )
00113 {
00114 typedef typename Polynomial::Scalar Scalar;
00115 typedef typename NumTraits<Scalar>::Real Real;
00116
00117 DenseIndex i=0;
00118 while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
00119 if( poly.size()-1 == i ){
00120 return Real(1); }
00121
00122 const Scalar inv_min_coeff = Scalar(1)/poly[i];
00123 Real cb(1);
00124 for( DenseIndex j=i+1; j<poly.size(); ++j ){
00125 cb += internal::abs(poly[j]*inv_min_coeff); }
00126 return Real(1)/cb;
00127 }
00128
00139 template <typename RootVector, typename Polynomial>
00140 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
00141 {
00142
00143 typedef typename Polynomial::Scalar Scalar;
00144
00145 poly.setZero( rv.size()+1 );
00146 poly[0] = -rv[0]; poly[1] = Scalar(1);
00147 for( DenseIndex i=1; i< rv.size(); ++i )
00148 {
00149 for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
00150 poly[0] = -rv[i]*poly[0];
00151 }
00152 }
00153
00154 }
00155
00156 #endif // EIGEN_POLYNOMIAL_UTILS_H