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_