core/vnl/vnl_polynomial.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_polynomial.h
00002 #ifndef vnl_polynomial_h_
00003 #define vnl_polynomial_h_
00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00005 #pragma interface
00006 #endif
00007 //:
00008 // \file
00009 // \brief Evaluation of univariate polynomials
00010 // Templated class (on the data type of the coefficients),
00011 // further very similar to the vnl_real_polynomial class,
00012 // except that it uses std::vector instead of vnl_vector as data container,
00013 // that the zero polynomial is represented by an empty vector,
00014 // and that the coefficients go in the other direction.
00015 //
00016 // Important note on the implementation choice (reversed coefficient vector
00017 // as opposed to the class vnl_real_npolynomial):
00018 // The choice made here is definitely the more natural one, since it makes
00019 // polynomials of different degrees much more naturally comparable, and hence
00020 // simplifies the implementation of e.g. operator+(). Indeed: even if the
00021 // degrees are different, the coefficients [i] of two polynomials are the ones
00022 // to be considered together since they both refer to $X^i$. Also, normalizing
00023 // the internal representation (in case the highest order coefficient is zero)
00024 // now just needs to pop_back() instead of shifting the coefficients vector.
00025 // In summary, the choice made here is both more natural and more performant.
00026 //
00027 // \author Peter Vanroose, ABIS Leuven.
00028 // \date  August 2011
00029 //
00030 // \verbatim
00031 //  Modifications
00032 //   20 Aug 2011 - Peter Vanroose - internal repr change: coeff vector reversed
00033 // \endverbatim
00034 
00035 #include <vcl_vector.h>
00036 #include <vcl_iosfwd.h>
00037 #include <vcl_cassert.h>
00038 
00039 //: Evaluation of polynomials.
00040 //  vnl_polynomial<T> represents a univariate polynomial with
00041 //  coefficients of datatype T, stored as a vector of values.
00042 //  This allows evaluation of the polynomial $p(X)$ at given values of $X$,
00043 //  and of its derivative $p'(X)$ or primitive function $\int p$.
00044 //
00045 //  The class also provides the common polynomial arithmetic, i.e.,
00046 //  + - *, and even long division (through operators / and %).
00047 //
00048 //  The coefficients (coeffs_) are stored as a vector, starting with
00049 //  the constant term. Hence coeffs_[n] is the coefficient of $X^n$,
00050 
00051 template <class T>
00052 class vnl_polynomial
00053 {
00054  public:
00055   //: Initialize the polynomial from its coefficients, lowest order first.
00056   // The polynomial is $ a[d] X^d + a[d-1] X^{d-1} + \cdots + a[0] = 0 $.
00057   // Note that this constructor expects the constant term coefficient first,
00058   // as opposed to the C array constructor!
00059   // An assertion makes sure that the input vector is in normalised form, i.e.,
00060   // that it is either empty or that the highest order coefficient is nonzero.
00061   vnl_polynomial(vcl_vector<T> const& a): coeffs_(a) { assert(a.begin()==a.end() || a.back() != T(0)); }
00062 
00063   //: Initialize polynomial from C array, highest order first.
00064   // The parameter \p len is the number of coefficients passed in,
00065   // which equals the degree plus one.
00066   // Note that this constructor expects the highest order coefficients first,
00067   // as opposed to the std::vector constructor!
00068   vnl_polynomial(T const* a, unsigned len) { assert(len==0 || *a != T(0)); while (len--) coeffs_.push_back(a[len]); }
00069 
00070   //: Initialize polynomial from single value, thus creating a monomial.
00071   // This is effectively an implicit cast from type T to class vnl_polynomial,
00072   // useful when adding or multiplying a polynomial with a number.
00073   vnl_polynomial(T const& a): coeffs_(1u, a) { if (a==T(0)) coeffs_.clear(); }
00074 
00075   //: Initialize polynomial of a given degree.
00076   // The default constructor initializes to the zero polynomial (which has degree -1).
00077   // but even with an explicit argument, the polynomial is the zero polynomial
00078   // (with non-compact storage) so it should always be used in conjunction with
00079   // operator[] for setting individual coefficients, at least coefficient [d].
00080   vnl_polynomial(int d=-1): coeffs_(d+1) { assert (d>=-1); }
00081 
00082   //: comparison operator
00083   bool operator==(vnl_polynomial<T> const& p) const { return p.coefficients() == coeffs_; }
00084 
00085   //: Returns negative of this polynomial
00086   vnl_polynomial<T> operator-() const;
00087 
00088   //: Returns polynomial which is sum of this with polynomial f
00089   vnl_polynomial<T> operator+(vnl_polynomial<T> const& f) const;
00090 
00091   //: Returns polynomial which is difference of this with polynomial f
00092   vnl_polynomial<T> operator-(vnl_polynomial<T> const& f) const { return operator+(-f); }
00093 
00094   //: Returns polynomial which is product of this with polynomial f
00095   vnl_polynomial<T> operator*(vnl_polynomial<T> const& f) const;
00096 
00097   //: Returns polynomial which is the result of the long division by polynomial f
00098   // Beware that this operation might not make sense for integral types T
00099   // if the highest order coefficient of f is not 1 or -1!
00100   vnl_polynomial<T> operator/(vnl_polynomial<T> const& f) const;
00101 
00102   //: Returns polynomial which is the remainder after a long division by polynomial f
00103   // Beware that this operation might not make sense for integral types T
00104   // if the highest order coefficient of f is not 1 or -1!
00105   vnl_polynomial<T> operator%(vnl_polynomial<T> const& f) const;
00106 
00107   vnl_polynomial<T>& operator+=(vnl_polynomial<T> const& f) { return *this = operator+(f); }
00108   vnl_polynomial<T>& operator-=(vnl_polynomial<T> const& f) { return *this = operator-(f); }
00109   vnl_polynomial<T>& operator*=(vnl_polynomial<T> const& f) { return *this = operator*(f); }
00110   vnl_polynomial<T>& operator/=(vnl_polynomial<T> const& f) { return *this = operator/(f); }
00111   vnl_polynomial<T>& operator%=(vnl_polynomial<T> const& f) { return *this = operator%(f); }
00112 
00113   //: Evaluate polynomial at value \p x
00114   T evaluate(T const& x) const;
00115 
00116   //: Return derivative of this polynomial
00117   vnl_polynomial<T> derivative() const;
00118 
00119   //: Evaluate derivative at value \p x
00120   T devaluate(T const& x) const { return derivative().evaluate(x); }
00121 
00122   //: Return primitive function (inverse derivative) of this polynomial
00123   // Since a primitive function is not unique, the one with constant term 0 is returned.
00124   // Beware that this operation might not make sense for integral types T!
00125   vnl_polynomial<T> primitive() const;
00126 
00127   //: Evaluate integral at \p x (assuming constant of integration is zero)
00128   // Beware that this operation might not make sense for integral types T!
00129   T evaluate_integral(T const& x) const { return primitive().evaluate(x); }
00130 
00131   //: Evaluate integral between \p x1 and \p x2
00132   // Beware that this operation might not make sense for integral types T!
00133   T evaluate_integral(T const& x1, T const& x2) const { return evaluate_integral(x2)-evaluate_integral(x1); }
00134 
00135   // Data Access---------------------------------------------------------------
00136 
00137   //: Return the degree (highest power of X) of the polynomial.
00138   // If the polynomial is zero, the degree is effectively -1.
00139   // Note that this method assumes a compactified representation, i.e., one
00140   // where the highest order coefficient is non-zero. Otherwise, the value
00141   // returned by degree() will be larger than the degree.
00142   int     degree() const { return int(coeffs_.size()) - 1; }
00143 
00144   //: Access to the polynomial coefficient of $X^i$
00145   T& operator [] (unsigned int i)       { assert(int(i)<=degree()); return coeffs_[i]; }
00146   //: Access to the polynomial coefficient of $X^i$
00147   T  operator [] (unsigned int i) const { assert(int(i)<=degree()); return coeffs_[i]; }
00148 
00149   //: Return the vector of coefficients
00150   const vcl_vector<T>& coefficients() const { return coeffs_; }
00151   //: Return the vector of coefficients
00152         vcl_vector<T>& coefficients()       { return coeffs_; }
00153 
00154   void set_coefficients(vcl_vector<T> const& coeffs) {coeffs_ = coeffs;}
00155 
00156   //: Print this polynomial to stream
00157   void print(vcl_ostream& os) const;
00158 
00159  protected:
00160   //: The coefficients of the polynomial.
00161   // coeffs_.front() is the const term.
00162   // coeffs_[n] is the coefficient of the $X^n$ term
00163   vcl_vector<T> coeffs_;
00164 };
00165 
00166 template <class T>
00167 vcl_ostream& operator<<(vcl_ostream& os, vnl_polynomial<T> const& p) { p.print(os); return os; }
00168 
00169 #define VNL_POLYNOMIAL_INSTANTIATE(T) extern "please #include vnl/vnl_polynomial.txx instead"
00170 
00171 #endif // vnl_polynomial_h_