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

details vigra/polynomial.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)