00001 // This is core/vnl/vnl_real_polynomial.h 00002 #ifndef vnl_real_polynomial_h_ 00003 #define vnl_real_polynomial_h_ 00004 #ifdef VCL_NEEDS_PRAGMA_INTERFACE 00005 #pragma interface 00006 #endif 00007 //: 00008 // \file 00009 // \brief Evaluation of real polynomials 00010 // \author Andrew W. Fitzgibbon, Oxford RRG 00011 // \date 06 Aug 96 00012 // 00013 // \verbatim 00014 // Modifications 00015 // 23 may 97, Peter Vanroose - "NO_COMPLEX" option added (until "complex" type is standardised) 00016 // 27/03/2001 Ian Scott and Tim Cootes - Added Binary IO 00017 // 27/03/2001 Ian Scott - Comments tidied up 00018 // 25/11/2001 Peter Vanroose - added operator==(), derivative(), primitive(), print() 00019 // 12/22/2004 Kongbin Kang - add structured comment for operator==() 00020 // \endverbatim 00021 00022 #include <vnl/vnl_vector.h> 00023 #include <vcl_complex.h> 00024 #include <vcl_iosfwd.h> 00025 #include <vcl_cassert.h> 00026 00027 //:Evaluation of real polynomials at real and complex points. 00028 // vnl_real_polynomial represents a univariate polynomial with real 00029 // coefficients, stored as a vector of doubles. This allows 00030 // evaluation of the polynomial $p(x)$ at given values of $x$, 00031 // or of its derivative $p'(x)$. 00032 // 00033 // The coefficients (coeffs_) are stored as a vnl_vector, where 00034 // coeffs_[n] is the coefficient of the x^(d-n) term, 00035 // where d is the degree of the polynomial. Otherwise said, 00036 // the coefficients are stored starting with the highest degree term. 00037 // 00038 // Roots may be extracted using the roots() method. 00039 class vnl_real_polynomial 00040 { 00041 public: 00042 //: Initialize polynomial. 00043 // The polynomial is $ a[0] x^d + a[1] x^{d-1} + \cdots + a[d] = 0 $. 00044 vnl_real_polynomial(vnl_vector<double> const & a): coeffs_(a) { 00045 if (a.empty()) { coeffs_.set_size(1); coeffs_(0)=0.0; } 00046 } 00047 00048 //: Initialize polynomial from C vector. 00049 // The parameter len is the number 00050 // of coefficients, one greater than the degree. 00051 vnl_real_polynomial(double const * a, unsigned len): coeffs_(a, len) { 00052 if (len==0) { coeffs_.set_size(1); coeffs_(0)=0.0; } 00053 } 00054 00055 //: Initialize polynomial from double. 00056 // Useful when adding or multiplying a polynomial and a number. 00057 vnl_real_polynomial(double a): coeffs_(1u, a) {} 00058 00059 //: Initialize polynomial of a given degree. 00060 vnl_real_polynomial(int d): coeffs_(d+1) { assert (d>=0); } 00061 00062 //: comparison operator 00063 bool operator==(vnl_real_polynomial const& p) const { return p.coefficients() == coeffs_; } 00064 00065 //: Evaluate polynomial at value x 00066 double evaluate(double x) const; 00067 00068 //: Evaluate integral at x (assuming constant of integration is zero) 00069 double evaluate_integral(double x) const; 00070 00071 //: Evaluate integral between x1 and x2 00072 double evaluate_integral(double x1, double x2) const; 00073 00074 //: Evaluate derivative at value x 00075 double devaluate(double x) const; 00076 00077 //: Evaluate polynomial at complex value x 00078 vcl_complex<double> evaluate(vcl_complex<double> const& x) const; 00079 00080 00081 //: Evaluate derivative at complex value x 00082 vcl_complex<double> devaluate(vcl_complex<double> const& x) const; 00083 00084 //: Return derivative of this polynomial 00085 vnl_real_polynomial derivative() const; 00086 00087 //: Return primitive function (inverse derivative) of this polynomial 00088 // Since a primitive function is not unique, the one with constant = 0 is returned 00089 vnl_real_polynomial primitive() const; 00090 00091 //: Add rhs to this and return *this 00092 vnl_real_polynomial& operator+=(vnl_real_polynomial const& rhs); 00093 00094 //: Subtract rhs from this and return *this 00095 vnl_real_polynomial& operator-=(vnl_real_polynomial const& rhs); 00096 00097 //: Multiply rhs with this and return *this 00098 vnl_real_polynomial& operator*=(vnl_real_polynomial const& rhs); 00099 00100 // Data Access--------------------------------------------------------------- 00101 00102 //: Return the degree (highest power of x) of the polynomial. 00103 int degree() const { return int(coeffs_.size()) - 1; } 00104 00105 //: Access to the polynomial coefficients 00106 double& operator [] (int i) { return coeffs_[i]; } 00107 //: Access to the polynomial coefficients 00108 double operator [] (int i) const { return coeffs_[i]; } 00109 00110 //: Return the vector of coefficients 00111 const vnl_vector<double>& coefficients() const { return coeffs_; } 00112 //: Return the vector of coefficients 00113 vnl_vector<double>& coefficients() { return coeffs_; } 00114 00115 void set_coefficients(vnl_vector<double> const& coeffs) {coeffs_ = coeffs;} 00116 00117 //: Print this polynomial to stream 00118 void print(vcl_ostream& os) const; 00119 00120 protected: 00121 //: The coefficients of the polynomial. 00122 // coeffs_.back() is the const term. 00123 // coeffs_[n] is the coefficient of the x^(d-n) term, 00124 // where d=coeffs_.size()-1 00125 // \invariant coeffs_size() >= 1; 00126 vnl_vector<double> coeffs_; 00127 }; 00128 00129 //: Returns polynomial which is sum of two polynomials f1(x)+f2(x) 00130 // \relatesalso vnl_real_polynomial 00131 vnl_real_polynomial operator+(const vnl_real_polynomial& f1, const vnl_real_polynomial& f2); 00132 00133 //: Returns polynomial which is different of two polynomials f1(x)-f2(x) 00134 // \relatesalso vnl_real_polynomial 00135 vnl_real_polynomial operator-(const vnl_real_polynomial& f1, const vnl_real_polynomial& f2); 00136 00137 //: Returns polynomial which is product of two polynomials f1(x)*f2(x) 00138 vnl_real_polynomial operator*(const vnl_real_polynomial& f1, const vnl_real_polynomial& f2); 00139 00140 //: Returns RMS difference between f1 and f2 over range [x1,x2] 00141 // $\frac1{\sqrt{|x_2-x_1|}}\,\sqrt{\int_{x_1}^{x_2}\left(f_1(x)-f_2(x)\right)^2\,dx}$ 00142 // \relatesalso vnl_real_polynomial 00143 double vnl_rms_difference(const vnl_real_polynomial& f1, const vnl_real_polynomial& f2, 00144 double x1, double x2); 00145 00146 #endif // vnl_real_polynomial_h_