Go to the documentation of this file.00001
00002 #ifndef vnl_polynomial_txx_
00003 #define vnl_polynomial_txx_
00004
00005 #include "vnl_polynomial.h"
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <vcl_iostream.h>
00015 #include <vcl_cassert.h>
00016
00017
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;
00025 for (++i; i!=coeffs_.end(); ++i) {
00026 acc += *i * xi;
00027 xi *= x;
00028 }
00029 return acc;
00030 }
00031
00032
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
00043 template <class T>
00044 vnl_polynomial<T> vnl_polynomial<T>::operator+(vnl_polynomial<T> const& f) const
00045 {
00046
00047 int d=degree(), d2=f.degree();
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
00052 while (sum.end() != sum.begin() && sum.back() == T(0)) sum.pop_back();
00053 return vnl_polynomial<T>(sum);
00054 }
00055
00056
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>();
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
00070
00071
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;
00076 assert (d2 >= 0 && f[d2] != T(0));
00077 if (d<0) return vnl_polynomial<T>();
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
00088
00089
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();
00097 vcl_vector<T> diff;
00098 for (int i=0; i<n; ++i) diff.push_back(coeffs_[i] - prod[i]);
00099
00100 while (diff.end() != diff.begin() && diff.back() == T(0)) diff.pop_back();
00101 return vnl_polynomial<T>(diff);
00102 }
00103
00104
00105 template <class T>
00106 vnl_polynomial<T> vnl_polynomial<T>::derivative() const
00107 {
00108 vcl_vector<T> cd;
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
00117
00118
00119 template <class T>
00120 vnl_polynomial<T> vnl_polynomial<T>::primitive() const
00121 {
00122 vcl_vector<T> cd;
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_