[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/polynomial.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2004 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00034 /*                                                                      */
00035 /************************************************************************/
00036 
00037 
00038 #ifndef VIGRA_POLYNOMIAL_HXX
00039 #define VIGRA_POLYNOMIAL_HXX
00040 
00041 #include <cmath>
00042 #include <complex>
00043 #include <algorithm>
00044 #include <iosfwd>
00045 #include "error.hxx"
00046 #include "mathutil.hxx"
00047 #include "numerictraits.hxx"
00048 #include "array_vector.hxx"
00049 
00050 namespace vigra {
00051 
00052 template <class T> class Polynomial;
00053 template <unsigned int MAXORDER, class T> class StaticPolynomial;
00054 
00055 /*****************************************************************/
00056 /*                                                               */
00057 /*                         PolynomialView                        */
00058 /*                                                               */
00059 /*****************************************************************/
00060 
00061 /** Polynomial interface for an externally managed array.
00062 
00063     The coefficient type <tt>T</tt> can be either a scalar or complex 
00064     (compatible to <tt>std::complex</tt>) type.
00065     
00066     \see vigra::Polynomial, vigra::StaticPolynomial, polynomialRoots()
00067 
00068     <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
00069     Namespace: vigra
00070     
00071     \ingroup Polynomials
00072 */
00073 template <class T>
00074 class PolynomialView
00075 {
00076   public:
00077     
00078         /** Coefficient type of the polynomial
00079         */
00080     typedef T         value_type;
00081 
00082         /** Promote type of <tt>value_type</tt>
00083             used for polynomial calculations
00084         */
00085     typedef typename NumericTraits<T>::RealPromote RealPromote;
00086 
00087         /** Scalar type associated with <tt>RealPromote</tt>
00088         */
00089     typedef typename NumericTraits<RealPromote>::ValueType Real;
00090 
00091         /** Complex type associated with <tt>RealPromote</tt>
00092         */
00093     typedef typename NumericTraits<RealPromote>::ComplexPromote Complex;
00094 
00095         /** Iterator for the coefficient sequence
00096         */
00097     typedef T *       iterator;
00098 
00099         /** Const iterator for the coefficient sequence
00100         */
00101     typedef T const * const_iterator;
00102     
00103     typedef Polynomial<Real> RealPolynomial;
00104     typedef Polynomial<Complex> ComplexPolynomial;
00105     
00106     
00107         /** Construct from a coefficient array of given <tt>order</tt>.
00108 
00109             The externally managed array must have length <tt>order+1</tt>
00110             and is interpreted as representing the polynomial:
00111             
00112             \code
00113             coeffs[0] + x * coeffs[1] + x * x * coeffs[2] + ...
00114             \endcode
00115             
00116             The coefficients are not copied, we only store a pointer to the 
00117             array.<tt>epsilon</tt> (default: 1.0e-14) determines the precision 
00118             of subsequent algorithms (especially root finding) performed on the
00119             polynomial.
00120         */
00121     PolynomialView(T * coeffs, unsigned int order, double epsilon = 1.0e-14)
00122     : coeffs_(coeffs),
00123       order_(order),
00124       epsilon_(epsilon)
00125     {}
00126     
00127         /// Access the coefficient of x^i
00128     T & operator[](unsigned int i)
00129         { return coeffs_[i]; }
00130     
00131         /// Access the coefficient of x^i
00132     T const & operator[](unsigned int i) const
00133         { return coeffs_[i]; }
00134     
00135         /** Evaluate the polynomial at the point <tt>v</tt> 
00136         
00137             Multiplication must be defined between the types
00138             <tt>V</tt> and <tt>PromoteTraits<T, V>::Promote</tt>.
00139             If both <tt>V</tt> and <tt>V</tt> are scalar, the result will
00140             be a scalar, otherwise it will be complex.
00141         */
00142     template <class V>
00143     typename PromoteTraits<T, V>::Promote
00144     operator()(V v) const;
00145     
00146         /** Differentiate the polynomial <tt>n</tt> times.
00147         */
00148     void differentiate(unsigned int n = 1);
00149     
00150         /** Deflate the polynomial at the root <tt>r</tt> with 
00151             the given <tt>multiplicity</tt>.
00152             
00153             The behavior of this function is undefined if <tt>r</tt>
00154             is not a root with at least the given multiplicity.
00155             This function calls forwardBackwardDeflate().
00156         */
00157     void deflate(T const & r, unsigned int multiplicity = 1);
00158     
00159         /** Forward-deflate the polynomial at the root <tt>r</tt>.
00160             
00161             The behavior of this function is undefined if <tt>r</tt>
00162             is not a root. Forward deflation is best if <tt>r</tt> is
00163             the biggest root (by magnitude).
00164         */
00165     void forwardDeflate(T const & v);
00166     
00167         /** Forward/backward eflate the polynomial at the root <tt>r</tt>.
00168             
00169             The behavior of this function is undefined if <tt>r</tt>
00170             is not a root. Combined forward/backward deflation is best 
00171             if <tt>r</tt> is an ontermediate root or we don't know
00172             <tt>r</tt>'s relation to the other roots of the polynomial.
00173         */
00174     void forwardBackwardDeflate(T v);
00175     
00176         /** Backward-deflate the polynomial at the root <tt>r</tt>.
00177             
00178             The behavior of this function is undefined if <tt>r</tt>
00179             is not a root. Backward deflation is best if <tt>r</tt> is
00180             the snallest root (by magnitude).
00181         */
00182     void backwardDeflate(T v);
00183     
00184         /** Deflate the polynomial with the complex conjugate roots 
00185             <tt>r</tt> and <tt>conj(r)</tt>.
00186             
00187             The behavior of this function is undefined if these are not
00188             roots.
00189         */
00190     void deflateConjugatePair(Complex const & v);
00191     
00192         /** Adjust the polynomial's order if the highest coefficients are near zero.
00193             The order is reduced as long as the absolute value does not exceed 
00194             the given \a epsilon.
00195         */
00196     void minimizeOrder(double epsilon = 0.0);    
00197     
00198         /** Normalize the polynomial, i.e. dived by the highest coefficient.
00199         */
00200     void normalize();
00201      
00202     void balance();
00203         
00204         /** Get iterator for the coefficient sequence.
00205         */
00206     iterator begin()
00207         { return coeffs_; }
00208     
00209         /** Get end iterator for the coefficient sequence.
00210         */
00211     iterator end()
00212         { return begin() + size(); }
00213     
00214         /** Get const_iterator for the coefficient sequence.
00215         */
00216     const_iterator begin() const
00217         { return coeffs_; }
00218     
00219         /** Get end const_iterator for the coefficient sequence.
00220         */
00221     const_iterator end() const
00222         { return begin() + size(); }
00223     
00224         /** Get length of the coefficient sequence (<tt>order() + 1</tt>).
00225         */
00226     unsigned int size() const
00227         { return order_ + 1; }
00228         
00229         /** Get order of the polynomial.
00230         */
00231     unsigned int order() const
00232         { return order_; }
00233         
00234         /** Get requested precision for polynomial algorithms 
00235             (especially root finding).
00236         */
00237     double epsilon() const
00238         { return epsilon_; }
00239         
00240         /** Set requested precision for polynomial algorithms 
00241             (especially root finding).
00242         */
00243     void setEpsilon(double eps)
00244         { epsilon_ = eps; }
00245 
00246   protected:
00247     PolynomialView(double epsilon = 1e-14)
00248     : coeffs_(0),
00249       order_(0),
00250       epsilon_(epsilon)
00251     {}
00252     
00253     void setCoeffs(T * coeffs, unsigned int order)
00254     {
00255         coeffs_ = coeffs;
00256         order_ = order;
00257     }
00258   
00259     T * coeffs_;
00260     unsigned int order_;
00261     double epsilon_;
00262 };
00263 
00264 template <class T>
00265 template <class U>
00266 typename PromoteTraits<T, U>::Promote
00267 PolynomialView<T>::operator()(U v) const
00268 {
00269     typename PromoteTraits<T, U>::Promote p(coeffs_[order_]);
00270     for(int i = order_ - 1; i >= 0; --i)
00271     {
00272        p = v * p + coeffs_[i];
00273     }
00274     return p;
00275 }
00276 
00277 /*
00278 template <class T>
00279 typename PolynomialView<T>::Complex 
00280 PolynomialView<T>::operator()(Complex const & v) const
00281 {
00282     Complex p(coeffs_[order_]);
00283     for(int i = order_ - 1; i >= 0; --i)
00284     {
00285        p = v * p + coeffs_[i];
00286     }
00287     return p;
00288 }
00289 */
00290 
00291 template <class T>
00292 void
00293 PolynomialView<T>::differentiate(unsigned int n)
00294 {
00295     if(n == 0)
00296         return;
00297     if(order_ == 0)
00298     {
00299         coeffs_[0] = 0.0;
00300         return;
00301     }
00302     for(unsigned int i = 1; i <= order_; ++i)
00303     {
00304         coeffs_[i-1] = double(i)*coeffs_[i];
00305     }
00306     --order_;
00307     if(n > 1)
00308         differentiate(n-1);
00309 }
00310 
00311 template <class T>
00312 void
00313 PolynomialView<T>::deflate(T const & v, unsigned int multiplicity)
00314 {
00315     vigra_precondition(order_ > 0,
00316         "PolynomialView<T>::deflate(): cannot deflate 0th order polynomial.");
00317     if(v == 0.0)
00318     {
00319         ++coeffs_;
00320         --order_;
00321     }
00322     else
00323     {
00324         // we use combined forward/backward deflation because
00325         // our initial guess seems to favour convergence to 
00326         // a root with magnitude near the median among all roots
00327         forwardBackwardDeflate(v);
00328     }
00329     if(multiplicity > 1)
00330         deflate(v, multiplicity-1);
00331 }
00332 
00333 template <class T>
00334 void
00335 PolynomialView<T>::forwardDeflate(T const & v)
00336 {
00337     for(int i = order_-1; i > 0; --i)
00338     {
00339         coeffs_[i] += v * coeffs_[i+1];
00340     }
00341     ++coeffs_;
00342     --order_;
00343 }
00344 
00345 template <class T>
00346 void
00347 PolynomialView<T>::forwardBackwardDeflate(T v)
00348 {
00349     unsigned int order2 = order_ / 2;
00350     T tmp = coeffs_[order_];
00351     for(unsigned int i = order_-1; i >= order2; --i)
00352     {
00353         T tmp1 = coeffs_[i];
00354         coeffs_[i] = tmp;
00355         tmp = tmp1 + v * tmp;
00356     }
00357     v = -1.0 / v;
00358     coeffs_[0] *= v;
00359     for(unsigned int i = 1; i < order2; ++i)
00360     {
00361         coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
00362     }
00363     --order_;
00364 }
00365 
00366 template <class T>
00367 void
00368 PolynomialView<T>::backwardDeflate(T v)
00369 {
00370     v = -1.0 / v;
00371     coeffs_[0] *= v;
00372     for(unsigned int i = 1; i < order_; ++i)
00373     {
00374         coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
00375     }
00376     --order_;
00377 }
00378 
00379 template <class T>
00380 void
00381 PolynomialView<T>::deflateConjugatePair(Complex const & v)
00382 {
00383     vigra_precondition(order_ > 1,
00384         "PolynomialView<T>::deflateConjugatePair(): cannot deflate 2 roots "
00385         "from 1st order polynomial.");
00386     Real a = 2.0*v.real();
00387     Real b = -sq(v.real()) - sq(v.imag());
00388     coeffs_[order_-1] += a * coeffs_[order_];
00389     for(int i = order_-2; i > 1; --i)
00390     {
00391         coeffs_[i] += a * coeffs_[i+1] + b*coeffs_[i+2];
00392     }
00393     coeffs_ += 2;
00394     order_ -= 2;
00395 }
00396     
00397 template <class T>
00398 void 
00399 PolynomialView<T>::minimizeOrder(double epsilon)
00400 {
00401     while(std::abs(coeffs_[order_]) <= epsilon && order_ > 0)
00402             --order_;
00403 }
00404 
00405 template <class T>
00406 void 
00407 PolynomialView<T>::normalize()
00408 {
00409     for(unsigned int i = 0; i<order_; ++i)
00410         coeffs_[i] /= coeffs_[order_];
00411     coeffs_[order_] = T(1.0);
00412 }
00413 
00414 template <class T>
00415 void 
00416 PolynomialView<T>::balance()
00417 {
00418     Real p0 = abs(coeffs_[0]), po = abs(coeffs_[order_]);
00419     Real norm = (p0 > 0.0)
00420                     ? VIGRA_CSTD::sqrt(p0*po) 
00421                     : po;
00422     for(unsigned int i = 0; i<=order_; ++i)
00423         coeffs_[i] /= norm;
00424 }
00425 
00426 /*****************************************************************/
00427 /*                                                               */
00428 /*                           Polynomial                          */
00429 /*                                                               */
00430 /*****************************************************************/
00431 
00432 /** Polynomial with internally managed array.
00433 
00434     Most interesting functionality is inherited from \ref vigra::PolynomialView.
00435 
00436     \see vigra::PolynomialView, vigra::StaticPolynomial, polynomialRoots()
00437 
00438     <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
00439     Namespace: vigra
00440     
00441     \ingroup Polynomials
00442 */
00443 template <class T>
00444 class Polynomial
00445 : public PolynomialView<T>
00446 {
00447     typedef PolynomialView<T> BaseType;
00448   public:
00449     typedef typename BaseType::Real    Real;
00450     typedef typename BaseType::Complex Complex;
00451     typedef Polynomial<Real>           RealPolynomial;
00452     typedef Polynomial<Complex>        ComplexPolynomial;
00453 
00454     typedef T         value_type;
00455     typedef T *       iterator;
00456     typedef T const * const_iterator;    
00457     
00458         /** Construct polynomial with given <tt>order</tt> and all coefficients
00459             set to zero (they can be set later using <tt>operator[]</tt>
00460             or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 
00461             the precision of subsequent algorithms (especially root finding) 
00462             performed on the polynomial.
00463         */
00464     Polynomial(unsigned int order = 0, double epsilon = 1.0e-14)
00465     : BaseType(epsilon),
00466       polynomial_(order + 1, T())
00467     {
00468         this->setCoeffs(&polynomial_[0], order);
00469     }
00470     
00471         /** Copy constructor
00472         */
00473     Polynomial(Polynomial const & p)
00474     : BaseType(p.epsilon()),
00475       polynomial_(p.begin(), p.end())
00476     {
00477         this->setCoeffs(&polynomial_[0], p.order());
00478     }
00479 
00480         /** Construct polynomial by copying the given coefficient sequence.
00481         */
00482     template <class ITER>
00483     Polynomial(ITER i, unsigned int order)
00484     : BaseType(),
00485       polynomial_(i, i + order + 1)
00486     {
00487         this->setCoeffs(&polynomial_[0], order);
00488     }
00489     
00490         /** Construct polynomial by copying the given coefficient sequence.
00491             Set <tt>epsilon</tt> (default: 1.0e-14) as 
00492             the precision of subsequent algorithms (especially root finding) 
00493             performed on the polynomial.
00494         */
00495     template <class ITER>
00496     Polynomial(ITER i, unsigned int order, double epsilon)
00497     : BaseType(epsilon),
00498       polynomial_(i, i + order + 1)
00499     {
00500         this->setCoeffs(&polynomial_[0], order);
00501     }
00502     
00503         /** Assigment
00504         */
00505     Polynomial & operator=(Polynomial const & p)
00506     {
00507         if(this == &p)
00508             return *this;
00509         ArrayVector<T> tmp(p.begin(), p.end());
00510         polynomial_.swap(tmp);
00511         this->setCoeffs(&polynomial_[0], p.order());
00512         this->epsilon_ = p.epsilon_;
00513         return *this;
00514     }
00515     
00516         /** Construct new polynomial representing the derivative of this
00517             polynomial.
00518         */
00519     Polynomial<T> getDerivative(unsigned int n = 1) const
00520     {
00521         Polynomial<T> res(*this);
00522         res.differentiate(n);
00523         return res;
00524     }
00525 
00526         /** Construct new polynomial representing this polynomial after
00527             deflation at the real root <tt>r</tt>.
00528         */
00529     Polynomial<T> 
00530     getDeflated(Real r) const
00531     {
00532         Polynomial<T> res(*this);
00533         res.deflate(r);
00534         return res;
00535     }
00536 
00537         /** Construct new polynomial representing this polynomial after
00538             deflation at the complex root <tt>r</tt>. The resulting
00539             polynomial will have complex coefficients, even if this
00540             polynomial had real ones.
00541         */
00542     Polynomial<Complex> 
00543     getDeflated(Complex const & r) const
00544     {
00545         Polynomial<Complex> res(this->begin(), this->order(), this->epsilon());
00546         res.deflate(r);
00547         return res;
00548     }
00549 
00550   protected:
00551     ArrayVector<T> polynomial_;
00552 };
00553 
00554 /*****************************************************************/
00555 /*                                                               */
00556 /*                        StaticPolynomial                       */
00557 /*                                                               */
00558 /*****************************************************************/
00559 
00560 /** Polynomial with internally managed array of static length.
00561 
00562     Most interesting functionality is inherited from \ref vigra::PolynomialView.
00563     This class differs from \ref vigra::Polynomial in that it allocates
00564     its memory statically which is much faster. Therefore, <tt>StaticPolynomial</tt>
00565     can only represent polynomials up to the given <tt>MAXORDER</tt>.
00566 
00567     \see vigra::PolynomialView, vigra::Polynomial, polynomialRoots()
00568 
00569     <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
00570     Namespace: vigra
00571     
00572     \ingroup Polynomials
00573 */
00574 template <unsigned int MAXORDER, class T>
00575 class StaticPolynomial
00576 : public PolynomialView<T>
00577 {
00578     typedef PolynomialView<T> BaseType;
00579     
00580   public:
00581     typedef typename BaseType::Real    Real;
00582     typedef typename BaseType::Complex Complex;
00583     typedef StaticPolynomial<MAXORDER, Real> RealPolynomial;
00584     typedef StaticPolynomial<MAXORDER, Complex> ComplexPolynomial;
00585 
00586     typedef T         value_type;
00587     typedef T *       iterator;
00588     typedef T const * const_iterator;
00589     
00590     
00591         /** Construct polynomial with given <tt>order <= MAXORDER</tt> and all 
00592             coefficients set to zero (they can be set later using <tt>operator[]</tt>
00593             or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 
00594             the precision of subsequent algorithms (especially root finding) 
00595             performed on the polynomial.
00596         */
00597     StaticPolynomial(unsigned int order = 0, double epsilon = 1.0e-14)
00598     : BaseType(epsilon)
00599     {
00600         vigra_precondition(order <= MAXORDER,
00601             "StaticPolynomial(): order exceeds MAXORDER.");
00602         std::fill_n(polynomial_, order+1, T());
00603         this->setCoeffs(polynomial_, order);
00604     }
00605     
00606         /** Copy constructor
00607         */
00608     StaticPolynomial(StaticPolynomial const & p)
00609     : BaseType(p.epsilon())
00610     {
00611         std::copy(p.begin(), p.end(), polynomial_);
00612         this->setCoeffs(polynomial_, p.order());
00613     }
00614 
00615         /** Construct polynomial by copying the given coefficient sequence.
00616             <tt>order <= MAXORDER</tt> is required.
00617         */
00618     template <class ITER>
00619     StaticPolynomial(ITER i, unsigned int order)
00620     : BaseType()
00621     {
00622         vigra_precondition(order <= MAXORDER,
00623             "StaticPolynomial(): order exceeds MAXORDER.");
00624         std::copy(i, i + order + 1, polynomial_);
00625         this->setCoeffs(polynomial_, order);
00626     }
00627     
00628         /** Construct polynomial by copying the given coefficient sequence.
00629             <tt>order <= MAXORDER</tt> is required. Set <tt>epsilon</tt> (default: 1.0e-14) as 
00630             the precision of subsequent algorithms (especially root finding) 
00631             performed on the polynomial.
00632         */
00633     template <class ITER>
00634     StaticPolynomial(ITER i, unsigned int order, double epsilon)
00635     : BaseType(epsilon)
00636     {
00637         vigra_precondition(order <= MAXORDER,
00638             "StaticPolynomial(): order exceeds MAXORDER.");
00639         std::copy(i, i + order + 1, polynomial_);
00640         this->setCoeffs(polynomial_, order);
00641     }
00642     
00643         /** Assigment.
00644         */
00645     StaticPolynomial & operator=(StaticPolynomial const & p)
00646     {
00647         if(this == &p)
00648             return *this;
00649         std::copy(p.begin(), p.end(), polynomial_);
00650         this->setCoeffs(polynomial_, p.order());
00651         this->epsilon_ = p.epsilon_;
00652         return *this;
00653     }
00654     
00655         /** Construct new polynomial representing the derivative of this
00656             polynomial.
00657         */
00658     StaticPolynomial getDerivative(unsigned int n = 1) const
00659     {
00660         StaticPolynomial res(*this);
00661         res.differentiate(n);
00662         return res;
00663     }
00664 
00665         /** Construct new polynomial representing this polynomial after
00666             deflation at the real root <tt>r</tt>.
00667         */
00668     StaticPolynomial 
00669     getDeflated(Real r) const
00670     {
00671         StaticPolynomial res(*this);
00672         res.deflate(r);
00673         return res;
00674     }
00675 
00676         /** Construct new polynomial representing this polynomial after
00677             deflation at the complex root <tt>r</tt>. The resulting
00678             polynomial will have complex coefficients, even if this
00679             polynomial had real ones.
00680         */
00681     StaticPolynomial<MAXORDER, Complex> 
00682     getDeflated(Complex const & r) const
00683     {
00684         StaticPolynomial<MAXORDER, Complex>  res(this->begin(), this->order(), this->epsilon());
00685         res.deflate(r);
00686         return res;
00687     }
00688     
00689     void setOrder(unsigned int order)
00690     {
00691         vigra_precondition(order <= MAXORDER,
00692             "taticPolynomial::setOrder(): order exceeds MAXORDER.");
00693         this->order_ = order;
00694     }
00695 
00696   protected:
00697     T polynomial_[MAXORDER+1];
00698 };
00699 
00700 /************************************************************/
00701 
00702 namespace detail {
00703 
00704 // replacement for complex division (some compilers have numerically
00705 // less stable implementations); code from python complexobject.c
00706 template <class T>
00707 std::complex<T> complexDiv(std::complex<T> const & a, std::complex<T> const & b)
00708 {
00709      const double abs_breal = b.real() < 0 ? -b.real() : b.real();
00710      const double abs_bimag = b.imag() < 0 ? -b.imag() : b.imag();
00711 
00712      if (abs_breal >= abs_bimag) 
00713      {
00714         /* divide tops and bottom by b.real() */
00715         if (abs_breal == 0.0) 
00716         {
00717             return std::complex<T>(a.real() / abs_breal, a.imag() / abs_breal);
00718         }
00719         else 
00720         {
00721             const double ratio = b.imag() / b.real();
00722             const double denom = b.real() + b.imag() * ratio;
00723             return std::complex<T>((a.real() + a.imag() * ratio) / denom,
00724                                    (a.imag() - a.real() * ratio) / denom);
00725         }
00726     }
00727     else 
00728     {
00729         /* divide tops and bottom by b.imag() */
00730         const double ratio = b.real() / b.imag();
00731         const double denom = b.real() * ratio + b.imag();
00732         return std::complex<T>((a.real() * ratio + a.imag()) / denom,
00733                                (a.imag() * ratio - a.real()) / denom);
00734     }
00735 }
00736 
00737 template <class T>
00738 std::complex<T> deleteBelowEpsilon(std::complex<T> const & x, double eps)
00739 {
00740     return std::abs(x.imag()) <= 2.0*eps*std::abs(x.real()) 
00741                 ? std::complex<T>(x.real())
00742                 : std::abs(x.real()) <= 2.0*eps*std::abs(x.imag()) 
00743                     ? std::complex<T>(NumericTraits<T>::zero(), x.imag())
00744                     :  x;
00745 }
00746 
00747 template <class POLYNOMIAL>
00748 typename POLYNOMIAL::value_type
00749 laguerreStartingGuess(POLYNOMIAL const & p)
00750 {
00751     double N = p.order();
00752     typename POLYNOMIAL::value_type centroid = -p[p.order()-1] / N / p[p.order()];
00753     double dist = VIGRA_CSTD::pow(std::abs(p(centroid) / p[p.order()]), 1.0 / N);
00754     return centroid + dist;
00755 }
00756 
00757 template <class POLYNOMIAL, class Complex>
00758 int laguerre1Root(POLYNOMIAL const & p, Complex & x, unsigned int multiplicity)
00759 {
00760     typedef typename NumericTraits<Complex>::ValueType Real;
00761     
00762     static double frac[] = {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
00763     int maxiter = 80, 
00764         count;
00765     double N = p.order();
00766     double eps  = p.epsilon(),
00767            eps2 = VIGRA_CSTD::sqrt(eps);
00768         
00769     if(multiplicity == 0)
00770         x = laguerreStartingGuess(p);
00771         
00772     bool mayTryDerivative = true;  // try derivative for multiple roots
00773     
00774     for(count = 0; count < maxiter; ++count)
00775     {
00776         // Horner's algorithm to calculate values of polynomial and its
00777         // first two derivatives and estimate error for current x
00778         Complex p0(p[p.order()]);
00779         Complex p1(0.0);
00780         Complex p2(0.0);
00781         Real ax    = std::abs(x);
00782         Real err = std::abs(p0);
00783         for(int i = p.order()-1; i >= 0; --i)
00784         {
00785             p2  = p2  * x  + p1;
00786             p1  = p1  * x  + p0;
00787             p0  = p0  * x  + p[i];
00788             err = err * ax + std::abs(p0);
00789         }
00790         p2 *= 2.0;
00791         err *= eps;
00792         Real ap0 = std::abs(p0);
00793         if(ap0 <= err)
00794         {
00795             break;  // converged
00796         }
00797         Complex g = complexDiv(p1, p0);
00798         Complex g2 = g * g;
00799         Complex h = g2 - complexDiv(p2, p0);
00800         // estimate root multiplicity according to Tien Chen
00801         if(g2 != 0.0)
00802         {
00803             multiplicity = (unsigned int)VIGRA_CSTD::floor(N / 
00804                                 (std::abs(N * complexDiv(h, g2) - 1.0) + 1.0) + 0.5);
00805             if(multiplicity < 1)
00806                 multiplicity = 1;
00807         }
00808         // improve accuracy of multiple roots on the derivative, as suggested by C. Bond
00809         // (do this only if we are already near the root, otherwise we may converge to 
00810         //  a different root of the derivative polynomial)
00811         if(mayTryDerivative && multiplicity > 1 && ap0 < eps2)
00812         {
00813             Complex x1 = x;
00814             int derivativeMultiplicity = laguerre1Root(p.getDerivative(), x1, multiplicity-1);
00815             if(derivativeMultiplicity && std::abs(p(x1)) < std::abs(p(x)))
00816             {
00817                 // successful search on derivative
00818                 x = x1;
00819                 return derivativeMultiplicity + 1;
00820             }
00821             else
00822             {
00823                 // unsuccessful search on derivative => don't do it again
00824                 mayTryDerivative = false;
00825             }
00826         }
00827         Complex sq = VIGRA_CSTD::sqrt((N - 1.0) * (N * h - g2));
00828         Complex gp = g + sq;
00829         Complex gm = g - sq;
00830         if(std::abs(gp) < std::abs(gm))
00831             gp = gm;
00832         Complex dx;
00833         if(gp != 0.0)
00834         {
00835             dx = complexDiv(Complex(N) , gp);
00836         }
00837         else
00838         {
00839             // re-initialisation trick due to Numerical Recipes
00840             dx = (1.0 + ax) * Complex(VIGRA_CSTD::cos(double(count)), VIGRA_CSTD::sin(double(count)));
00841         }
00842         Complex x1 = x - dx;
00843 
00844         if(x1 - x == 0.0)
00845         {
00846             break;  // converged
00847         }
00848         if((count + 1) % 10)
00849             x = x1;
00850         else
00851             // cycle breaking trick according to Numerical Recipes
00852             x = x - frac[(count+1)/10] * dx;
00853     }
00854     return count < maxiter ? 
00855         multiplicity : 
00856         0;
00857 }
00858 
00859 template <class Real>
00860 struct PolynomialRootCompare
00861 {
00862     Real epsilon;
00863 
00864     PolynomialRootCompare(Real eps)
00865     : epsilon(eps)
00866     {}
00867     
00868     template <class T>
00869     bool operator()(T const & l, T const & r)
00870     {
00871         return closeAtTolerance(l.real(), r.real(), epsilon)
00872                      ? l.imag() < r.imag()
00873                      : l.real() < r.real();
00874     }
00875 };
00876 
00877 } // namespace detail 
00878 
00879 /** \addtogroup Polynomials Polynomials and root determination
00880 
00881     Classes to represent polynomials and functions to find polynomial roots.
00882 */
00883 //@{
00884 
00885 /*****************************************************************/
00886 /*                                                               */
00887 /*                         polynomialRoots                       */
00888 /*                                                               */
00889 /*****************************************************************/
00890 
00891 /** Determine the roots of the polynomial <tt>poriginal</tt>.
00892 
00893     The roots are appended to the vector <tt>roots</tt>, with optional root
00894     polishing as specified by <tt>polishRoots</tt> (default: do polishing). The function uses an 
00895     improved version of Laguerre's algorithm. The improvements are as follows:
00896     
00897     <ul>
00898     <li>It uses a clever initial guess for the iteration, according to a proposal by Tien Chen</li>
00899     <li>It estimates each root's multiplicity, again according to Tien Chen, and reduces multiplicity
00900         by switching to the polynomial's derivative (which has the same root, with multiplicity
00901         reduced by one), as proposed by C. Bond.</li>
00902     </ul>
00903     
00904     The algorithm has been successfully used for polynomials up to order 80.
00905     The function stops and returns <tt>false</tt> if an iteration fails to converge within 
00906     80 steps. The type <tt>POLYNOMIAL</tt> must be compatible to 
00907     \ref vigra::PolynomialView, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
00908     with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>.
00909 
00910     <b> Declaration:</b>
00911 
00912     pass arguments explicitly:
00913     \code
00914     namespace vigra {
00915         template <class POLYNOMIAL, class VECTOR>
00916         bool 
00917         polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots = true);
00918     }
00919     \endcode
00920 
00921 
00922     <b> Usage:</b>
00923 
00924         <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
00925         Namespace: vigra
00926 
00927     \code
00928     // encode the polynomial  x^4 - 1
00929     Polynomial<double> poly(4);
00930     poly[0] = -1.0;
00931     poly[4] =  1.0;
00932 
00933     ArrayVector<std::complex<double> > roots;
00934     polynomialRoots(poly, roots);
00935     \endcode
00936 
00937     \see polynomialRootsEigenvalueMethod()
00938 */
00939 template <class POLYNOMIAL, class VECTOR>
00940 bool polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots)
00941 {
00942     typedef typename POLYNOMIAL::value_type T;
00943     typedef typename POLYNOMIAL::Real    Real;
00944     typedef typename POLYNOMIAL::Complex Complex;
00945     typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial;
00946     
00947     double eps  = poriginal.epsilon();
00948 
00949     WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
00950     p.minimizeOrder();
00951     if(p.order() == 0)
00952         return true;
00953 
00954     Complex x = detail::laguerreStartingGuess(p);
00955     
00956     unsigned int multiplicity = 1;
00957     bool triedConjugate = false;
00958     
00959     // handle the high order cases
00960     while(p.order() > 2)
00961     {
00962         p.balance();
00963         
00964         // find root estimate using Laguerre's method on deflated polynomial p;
00965         // zero return indicates failure to converge
00966         multiplicity = detail::laguerre1Root(p, x, multiplicity);
00967     
00968         if(multiplicity == 0)
00969             return false;
00970         // polish root on original polynomial poriginal;
00971         // zero return indicates failure to converge
00972         if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
00973             return false;
00974         x = detail::deleteBelowEpsilon(x, eps);
00975         roots.push_back(x);
00976         p.deflate(x);
00977         // determine the next starting guess
00978         if(multiplicity > 1)
00979         {
00980             // probably multiple root => keep current root as starting guess
00981             --multiplicity;
00982             triedConjugate = false;
00983         }
00984         else
00985         {
00986             // need a new starting guess
00987             if(x.imag() != 0.0 && !triedConjugate)
00988             {
00989                 // if the root is complex and we don't already have 
00990                 // the conjugate root => try the conjugate as starting guess
00991                 triedConjugate = true;
00992                 x = conj(x);
00993             }
00994             else
00995             {
00996                 // otherwise generate new starting guess
00997                 triedConjugate = false;
00998                 x = detail::laguerreStartingGuess(p);
00999             }
01000         }
01001     }
01002     
01003     // handle the low order cases
01004     if(p.order() == 2)
01005     {
01006         Complex a = p[2];
01007         Complex b = p[1];
01008         Complex c = p[0];
01009         Complex b2 = std::sqrt(b*b - 4.0*a*c);
01010         Complex q;
01011         if((conj(b)*b2).real() >= 0.0)
01012             q = -0.5 * (b + b2);
01013         else
01014             q = -0.5 * (b - b2);
01015         x = detail::complexDiv(q, a);
01016         if(polishRoots)
01017             detail::laguerre1Root(poriginal, x, 1);
01018         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01019         x = detail::complexDiv(c, q);
01020         if(polishRoots)
01021             detail::laguerre1Root(poriginal, x, 1);
01022         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01023     }
01024     else if(p.order() == 1)
01025     {
01026         x = detail::complexDiv(-p[0], p[1]);
01027         if(polishRoots)
01028             detail::laguerre1Root(poriginal, x, 1);
01029         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01030     }
01031     std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
01032     return true;
01033 }
01034 
01035 template <class POLYNOMIAL, class VECTOR>
01036 inline bool 
01037 polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
01038 {
01039     return polynomialRoots(poriginal, roots, true);
01040 }
01041 
01042 /** Determine the real roots of the polynomial <tt>p</tt>.
01043 
01044     This function simply calls \ref polynomialRoots() and than throws away all complex roots.
01045     Accordingly, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
01046     with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>.
01047 
01048     <b> Declaration:</b>
01049 
01050     pass arguments explicitly:
01051     \code
01052     namespace vigra {
01053         template <class POLYNOMIAL, class VECTOR>
01054         bool 
01055         polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots = true);
01056     }
01057     \endcode
01058 
01059 
01060     <b> Usage:</b>
01061 
01062         <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
01063         Namespace: vigra
01064 
01065     \code
01066     // encode the polynomial  x^4 - 1
01067     Polynomial<double> poly(4);
01068     poly[0] = -1.0;
01069     poly[4] =  1.0;
01070 
01071     ArrayVector<double> roots;
01072     polynomialRealRoots(poly, roots);
01073     \endcode
01074 
01075     \see polynomialRealRootsEigenvalueMethod()
01076 */
01077 template <class POLYNOMIAL, class VECTOR>
01078 bool polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots)
01079 {
01080     typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
01081     ArrayVector<Complex> croots;
01082     if(!polynomialRoots(p, croots, polishRoots))
01083         return false;
01084     for(unsigned int i = 0; i < croots.size(); ++i)
01085         if(croots[i].imag() == 0.0)
01086             roots.push_back(croots[i].real());
01087     return true;
01088 }
01089 
01090 template <class POLYNOMIAL, class VECTOR>
01091 inline bool 
01092 polynomialRealRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
01093 {
01094     return polynomialRealRoots(poriginal, roots, true);
01095 }
01096 
01097 //@}
01098 
01099 } // namespace vigra
01100 
01101 namespace std {
01102 
01103 template <class T>
01104 ostream & operator<<(ostream & o, vigra::PolynomialView<T> const & p)
01105 {
01106     for(unsigned int k=0; k < p.order(); ++k)
01107         o << p[k] << " ";
01108     o << p[p.order()];
01109     return o;
01110 }
01111 
01112 } // namespace std 
01113 
01114 #endif // VIGRA_POLYNOMIAL_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)