src/Polynomials/PolynomialUtils.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_POLYNOMIAL_UTILS_H
0011 #define EIGEN_POLYNOMIAL_UTILS_H
0012
0013 namespace Eigen {
0014
0015 /** \ingroup Polynomials_Module
0016  * \returns the evaluation of the polynomial at x using Horner algorithm.
0017  *
0018  * \param[in] poly : the vector of coefficients of the polynomial ordered
0019  *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
0020  *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
0021  * \param[in] x : the value to evaluate the polynomial at.
0022  *
0023  * \note for stability:
0024  *   \f$ |x| \le 1 \f$
0025  */
0026 template <typename Polynomials, typename T>
0027 inline
0028 T poly_eval_horner( const Polynomials& poly, const T& x )
0029 {
0030   T val=poly[poly.size()-1];
0031   for(DenseIndex i=poly.size()-2; i>=0; --i ){
0032     val = val*x + poly[i]; }
0033   return val;
0034 }
0035
0036 /** \ingroup Polynomials_Module
0037  * \returns the evaluation of the polynomial at x using stabilized Horner algorithm.
0038  *
0039  * \param[in] poly : the vector of coefficients of the polynomial ordered
0040  *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
0041  *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
0042  * \param[in] x : the value to evaluate the polynomial at.
0043  */
0044 template <typename Polynomials, typename T>
0045 inline
0046 T poly_eval( const Polynomials& poly, const T& x )
0047 {
0048   typedef typename NumTraits<T>::Real Real;
0049
0050   if( numext::abs2( x ) <= Real(1) ){
0051     return poly_eval_horner( poly, x ); }
0052   else
0053   {
0054     T val=poly[0];
0055     T inv_x = T(1)/x;
0056     for( DenseIndex i=1; i<poly.size(); ++i ){
0057       val = val*inv_x + poly[i]; }
0058
0059     return numext::pow(x,(T)(poly.size()-1)) * val;
0060   }
0061 }
0062
0063 /** \ingroup Polynomials_Module
0064  * \returns a maximum bound for the absolute value of any root of the polynomial.
0065  *
0066  * \param[in] poly : the vector of coefficients of the polynomial ordered
0067  *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
0068  *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
0069  *
0070  *  \pre
0071  *   the leading coefficient of the input polynomial poly must be non zero
0072  */
0073 template <typename Polynomial>
0074 inline
0075 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_max_bound( const Polynomial& poly )
0076 {
0077   using std::abs;
0078   typedef typename Polynomial::Scalar Scalar;
0079   typedef typename NumTraits<Scalar>::Real Real;
0080
0081   eigen_assert( Scalar(0) != poly[poly.size()-1] );
0082   const Scalar inv_leading_coeff = Scalar(1)/poly[poly.size()-1];
0083   Real cb(0);
0084
0085   for( DenseIndex i=0; i<poly.size()-1; ++i ){
0086     cb += abs(poly[i]*inv_leading_coeff); }
0087   return cb + Real(1);
0088 }
0089
0090 /** \ingroup Polynomials_Module
0091  * \returns a minimum bound for the absolute value of any non zero root of the polynomial.
0092  * \param[in] poly : the vector of coefficients of the polynomial ordered
0093  *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
0094  *  e.g. \f$ 1 + 3x^2 \f$ is stored as a vector \f$ [ 1, 0, 3 ] \f$.
0095  */
0096 template <typename Polynomial>
0097 inline
0098 typename NumTraits<typename Polynomial::Scalar>::Real cauchy_min_bound( const Polynomial& poly )
0099 {
0100   using std::abs;
0101   typedef typename Polynomial::Scalar Scalar;
0102   typedef typename NumTraits<Scalar>::Real Real;
0103
0104   DenseIndex i=0;
0105   while( i<poly.size()-1 && Scalar(0) == poly(i) ){ ++i; }
0106   if( poly.size()-1 == i ){
0107     return Real(1); }
0108
0109   const Scalar inv_min_coeff = Scalar(1)/poly[i];
0110   Real cb(1);
0111   for( DenseIndex j=i+1; j<poly.size(); ++j ){
0112     cb += abs(poly[j]*inv_min_coeff); }
0113   return Real(1)/cb;
0114 }
0115
0116 /** \ingroup Polynomials_Module
0117  * Given the roots of a polynomial compute the coefficients in the
0118  * monomial basis of the monic polynomial with same roots and minimal degree.
0119  * If RootVector is a vector of complexes, Polynomial should also be a vector
0120  * of complexes.
0121  * \param[in] rv : a vector containing the roots of a polynomial.
0122  * \param[out] poly : the vector of coefficients of the polynomial ordered
0123  *  by degrees i.e. poly[i] is the coefficient of degree i of the polynomial
0124  *  e.g. \f$ 3 + x^2 \f$ is stored as a vector \f$ [ 3, 0, 1 ] \f$.
0125  */
0126 template <typename RootVector, typename Polynomial>
0127 void roots_to_monicPolynomial( const RootVector& rv, Polynomial& poly )
0128 {
0129
0130   typedef typename Polynomial::Scalar Scalar;
0131
0132   poly.setZero( rv.size()+1 );
0133   poly[0] = -rv[0]; poly[1] = Scalar(1);
0134   for( DenseIndex i=1; i< rv.size(); ++i )
0135   {
0136     for( DenseIndex j=i+1; j>0; --j ){ poly[j] = poly[j-1] - rv[i]*poly[j]; }
0137     poly[0] = -rv[i]*poly[0];
0138   }
0139 }
0140
0141 } // end namespace Eigen
0142
0143 #endif // EIGEN_POLYNOMIAL_UTILS_H