blob: 6af8613b6c1ae5f6b8c38b74e5eb99e6388f0af5 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2010 Manuel Yguel // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_POLYNOMIAL_UTILS_H #define EIGEN_POLYNOMIAL_UTILS_H #include "./InternalHeaderCheck.h" namespace Eigen { /** \ingroup Polynomials_Module * \returns the evaluation of the polynomial at x using Horner algorithm. * * \param[in] poly : the vector of coefficients of the polynomial ordered * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial * e.g. \f\$ 1 + 3x^2 \f\$ is stored as a vector \f\$ [ 1, 0, 3 ] \f\$. * \param[in] x : the value to evaluate the polynomial at. * * \note for stability: * \f\$ |x| \le 1 \f\$ */ template inline T poly_eval_horner( const Polynomials& poly, const T& x ) { T val=poly[poly.size()-1]; for(DenseIndex i=poly.size()-2; i>=0; --i ){ val = val*x + poly[i]; } return val; } /** \ingroup Polynomials_Module * \returns the evaluation of the polynomial at x using stabilized Horner algorithm. * * \param[in] poly : the vector of coefficients of the polynomial ordered * by degrees i.e. poly[i] is the coefficient of degree i of the polynomial * e.g. \f\$ 1 + 3x^2 \f\$ is stored as a vector \f\$ [ 1, 0, 3 ] \f\$. * \param[in] x : the value to evaluate the polynomial at. */ template inline T poly_eval( const Polynomials& poly, const T& x ) { typedef typename NumTraits::Real Real; if( numext::abs2( x ) <= Real(1) ){ return poly_eval_horner( poly, x ); } else { T val=poly[0]; T inv_x = T(1)/x; for( DenseIndex i=1; i inline typename NumTraits::Real cauchy_max_bound( const Polynomial& poly ) { using std::abs; typedef typename Polynomial::Scalar Scalar; typedef typename NumTraits::Real Real; eigen_assert( Scalar(0) != poly[poly.size()-1] ); const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1]; Real cb(0); for( DenseIndex i=0; i inline typename NumTraits::Real cauchy_min_bound( const Polynomial& poly ) { using std::abs; typedef typename Polynomial::Scalar Scalar; typedef typename NumTraits::Real Real; DenseIndex i=0; while( i void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly ) { typedef typename Polynomial::Scalar Scalar; poly.setZero( rv.size()+1 ); poly[0] = -rv[0]; poly[1] = Scalar(1); for( DenseIndex i=1; i< rv.size(); ++i ) { for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; } poly[0] = -rv[i]*poly[0]; } } } // end namespace Eigen #endif // EIGEN_POLYNOMIAL_UTILS_H