core/vnl/vnl_polynomial.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_polynomial.txx
00002 #ifndef vnl_polynomial_txx_
00003 #define vnl_polynomial_txx_
00004 
00005 #include "vnl_polynomial.h"
00006 //:
00007 // \file
00008 // \brief Evaluation of polynomials - the implementation
00009 // Templated class (on the data type of the coefficients),
00010 // further very similar to the vnl_real_polynomial class.
00011 // \author Peter Vanroose, ABIS Leuven.
00012 // \date  August 2011
00013 
00014 #include <vcl_iostream.h>
00015 #include <vcl_cassert.h>
00016 
00017 //: Evaluate polynomial at value x
00018 template <class T>
00019 T vnl_polynomial<T>::evaluate(T const& x) const
00020 {
00021   typename vcl_vector<T>::const_iterator i = coeffs_.begin();
00022   if (i == coeffs_.end()) return T(0);
00023   T acc = *i;
00024   T xi = x; // will be x^i
00025   for (++i; i!=coeffs_.end(); ++i) {
00026     acc += *i * xi;
00027     xi *= x;
00028   }
00029   return acc;
00030 }
00031 
00032 //: Returns negative of this polynomial
00033 template <class T>
00034 vnl_polynomial<T> vnl_polynomial<T>::operator-() const
00035 {
00036   vcl_vector<T> neg = coeffs_;
00037   typename vcl_vector<T>::iterator i = neg.begin();
00038   for (; i!=neg.end(); ++i) *i = - *i;
00039   return vnl_polynomial<T>(neg);
00040 }
00041 
00042 //: Returns polynomial which is sum of this with polynomial f
00043 template <class T>
00044 vnl_polynomial<T> vnl_polynomial<T>::operator+(vnl_polynomial<T> const& f) const
00045 {
00046   // Degree of result is (at most) the maximum of the two input degrees:
00047   int d=degree(), d2=f.degree(); // any or both of these might be -1 !
00048   vcl_vector<T> sum = coeffs_;
00049   for (int i=0;i<=d&&i<=d2;++i) sum[i]+=f[i];
00050   for (int i=d+1;i<=d2;++i) sum.push_back(f[i]);
00051   // normalise the result, viz. such that the highest order coefficient is zero:
00052   while (sum.end() != sum.begin() && sum.back() == T(0)) sum.pop_back();
00053   return vnl_polynomial<T>(sum);
00054 }
00055 
00056 //: Returns polynomial which is product of this with polynomial f
00057 template <class T>
00058 vnl_polynomial<T> vnl_polynomial<T>::operator*(vnl_polynomial<T> const& f) const
00059 {
00060   int d1=degree(), d2=f.degree(), d = d1+d2;
00061   if (d1<0 || d2<0) return vnl_polynomial<T>(); // one of the factors is zero
00062   vcl_vector<T> prod(d+1, T(0));
00063   for (int i=0;i<=d1;++i)
00064     for (int j=0;j<=d2;++j)
00065       prod[i+j] += coeffs_[i]*f[j];
00066   return vnl_polynomial<T>(prod);
00067 }
00068 
00069 //: Returns polynomial which is the result of a long division by f
00070 // Beware that this operation might not make sense for integral types T
00071 // if the highest order coefficient of f is not 1 or -1!
00072 template <class T>
00073 vnl_polynomial<T> vnl_polynomial<T>::operator/(vnl_polynomial<T> const& f) const
00074 {
00075   int d1=degree(), d2=f.degree(), d=d1-d2; // d will be the degree of the quotient
00076   assert (d2 >= 0 && f[d2] != T(0)); // denominator should not be zero
00077   if (d<0) return vnl_polynomial<T>(); // nominator is zero, or denominator has higher degree than nominator
00078   vcl_vector<T> quot;
00079   for (int i=0;i<=d;++i) {
00080     T acc = coeffs_[d1-i];
00081     for (int j=0;j<d2&&j<i;++j) acc -= quot[j] * f[d2-j-1];
00082     quot.insert(quot.begin(), 1, acc/f[d2]);
00083   }
00084   return vnl_polynomial<T>(quot);
00085 }
00086 
00087 //: Returns polynomial which is the remainder after long division by f
00088 // Beware that this operation might not make sense for integral types T
00089 // if the highest order coefficient of f is not 1 or -1!
00090 template <class T>
00091 vnl_polynomial<T> vnl_polynomial<T>::operator%(vnl_polynomial<T> const& f) const
00092 {
00093   vnl_polynomial<T> quot = operator/(f);
00094   if (quot.degree() < 0) return *this;
00095   vnl_polynomial<T> prod = f * quot;
00096   int n=f.degree(); // size of the result, i.e., one more than degree of the result
00097   vcl_vector<T> diff;
00098   for (int i=0; i<n; ++i) diff.push_back(coeffs_[i] - prod[i]);
00099   // normalise the result, viz. such that the highest order coefficient is zero:
00100   while (diff.end() != diff.begin() && diff.back() == T(0)) diff.pop_back();
00101   return vnl_polynomial<T>(diff);
00102 }
00103 
00104 //: Return derivative of this polynomial
00105 template <class T>
00106 vnl_polynomial<T> vnl_polynomial<T>::derivative() const
00107 {
00108   vcl_vector<T> cd; // will be one shorter than coeffs_
00109   typename vcl_vector<T>::const_iterator i = coeffs_.begin();
00110   T n = T(1);
00111   for (++i; i!=coeffs_.end(); ++i,++n)
00112     cd.push_back(*i * n);
00113   return vnl_polynomial<T>(cd);
00114 }
00115 
00116 //: Return primitive function (inverse derivative) of this polynomial
00117 // Since a primitive function is not unique, the one with constant = 0 is returned
00118 // Beware that this operation might not make sense for integral types T!
00119 template <class T>
00120 vnl_polynomial<T> vnl_polynomial<T>::primitive() const
00121 {
00122   vcl_vector<T> cd; // will be one longer than coeffs_
00123   T n = T(0);
00124   cd.push_back(n);
00125   typename vcl_vector<T>::const_iterator i = coeffs_.begin();
00126   for (++n; i!=coeffs_.end(); ++i,++n)
00127     cd.push_back(*i / n);
00128   return vnl_polynomial<T>(cd);
00129 }
00130 
00131 template <class T>
00132 void vnl_polynomial<T>::print(vcl_ostream& os) const
00133 {
00134   bool first_coeff = true;
00135 
00136   for (int i=degree(); i >= 0; --i) {
00137     if (coeffs_[i] == T(0)) continue;
00138     os << ' ';
00139     if (coeffs_[i] > T(0) && !first_coeff) os << '+';
00140     if (i==0)                     os << coeffs_[i];
00141     else if (coeffs_[i] == -T(1)) os << '-';
00142     else if (coeffs_[i] != T(1))  os << coeffs_[i] << ' ';
00143     if (i == 1)                   os << 'X';
00144     else if (i != 0)              os << "X^" << i;
00145     first_coeff = false;
00146   }
00147   if (first_coeff) os << " 0";
00148 }
00149 
00150 #undef VNL_POLYNOMIAL_INSTANTIATE
00151 #define VNL_POLYNOMIAL_INSTANTIATE(T) \
00152 template class vnl_polynomial<T >; \
00153 template vcl_ostream& operator<<(vcl_ostream& os, vnl_polynomial<T > const&)
00154 
00155 #endif // vnl_polynomial_txx_