core/vnl/vnl_real_polynomial.h
Go to the documentation of this file.
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_