[ 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.3, Aug 18 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     void balance();
00188         
00189         /** Get iterator for the coefficient sequence.
00190         */
00191     iterator begin()
00192         { return coeffs_; }
00193     
00194         /** Get end iterator for the coefficient sequence.
00195         */
00196     iterator end()
00197         { return begin() + size(); }
00198     
00199         /** Get const_iterator for the coefficient sequence.
00200         */
00201     const_iterator begin() const
00202         { return coeffs_; }
00203     
00204         /** Get end const_iterator for the coefficient sequence.
00205         */
00206     const_iterator end() const
00207         { return begin() + size(); }
00208     
00209         /** Get length of the coefficient sequence (<tt>order() + 1</tt>).
00210         */
00211     unsigned int size() const
00212         { return order_ + 1; }
00213         
00214         /** Get order of the polynomial.
00215         */
00216     unsigned int order() const
00217         { return order_; }
00218         
00219         /** Get requested precision for polynomial algorithms 
00220             (especially root finding).
00221         */
00222     double epsilon() const
00223         { return epsilon_; }
00224         
00225         /** Set requested precision for polynomial algorithms 
00226             (especially root finding).
00227         */
00228     void setEpsilon(double eps)
00229         { epsilon_ = eps; }
00230 
00231   protected:
00232     PolynomialView(double epsilon = 1e-14)
00233     : coeffs_(0),
00234       order_(0),
00235       epsilon_(epsilon)
00236     {}
00237     
00238     void setCoeffs(T * coeffs, unsigned int order)
00239     {
00240         coeffs_ = coeffs;
00241         order_ = order;
00242     }
00243   
00244     T * coeffs_;
00245     unsigned int order_;
00246     double epsilon_;
00247 };
00248 
00249 template <class T>
00250 template <class U>
00251 typename PromoteTraits<T, U>::Promote
00252 PolynomialView<T>::operator()(U v) const
00253 {
00254     typename PromoteTraits<T, U>::Promote p(coeffs_[order_]);
00255     for(int i = order_ - 1; i >= 0; --i)
00256     {
00257        p = v * p + coeffs_[i];
00258     }
00259     return p;
00260 }
00261 
00262 /*
00263 template <class T>
00264 typename PolynomialView<T>::Complex 
00265 PolynomialView<T>::operator()(Complex const & v) const
00266 {
00267     Complex p(coeffs_[order_]);
00268     for(int i = order_ - 1; i >= 0; --i)
00269     {
00270        p = v * p + coeffs_[i];
00271     }
00272     return p;
00273 }
00274 */
00275 
00276 template <class T>
00277 void
00278 PolynomialView<T>::differentiate(unsigned int n)
00279 {
00280     if(n == 0)
00281         return;
00282     if(order_ == 0)
00283     {
00284         coeffs_[0] = 0.0;
00285         return;
00286     }
00287     for(unsigned int i = 1; i <= order_; ++i)
00288     {
00289         coeffs_[i-1] = double(i)*coeffs_[i];
00290     }
00291     --order_;
00292     if(n > 1)
00293         differentiate(n-1);
00294 }
00295 
00296 template <class T>
00297 void
00298 PolynomialView<T>::deflate(T const & v, unsigned int multiplicity)
00299 {
00300     vigra_precondition(order_ > 0,
00301         "PolynomialView<T>::deflate(): cannot deflate 0th order polynomial.");
00302     if(v == 0.0)
00303     {
00304         ++coeffs_;
00305         --order_;
00306     }
00307     else
00308     {
00309         // we use combined forward/backward deflation because
00310         // our initial guess seems to favour convergence to 
00311         // a root with magnitude near the median among all roots
00312         forwardBackwardDeflate(v);
00313     }
00314     if(multiplicity > 1)
00315         deflate(v, multiplicity-1);
00316 }
00317 
00318 template <class T>
00319 void
00320 PolynomialView<T>::forwardDeflate(T const & v)
00321 {
00322     for(int i = order_-1; i > 0; --i)
00323     {
00324         coeffs_[i] += v * coeffs_[i+1];
00325     }
00326     ++coeffs_;
00327     --order_;
00328 }
00329 
00330 template <class T>
00331 void
00332 PolynomialView<T>::forwardBackwardDeflate(T v)
00333 {
00334     unsigned int order2 = order_ / 2;
00335     T tmp = coeffs_[order_];
00336     for(unsigned int i = order_-1; i >= order2; --i)
00337     {
00338         T tmp1 = coeffs_[i];
00339         coeffs_[i] = tmp;
00340         tmp = tmp1 + v * tmp;
00341     }
00342     v = -1.0 / v;
00343     coeffs_[0] *= v;
00344     for(unsigned int i = 1; i < order2; ++i)
00345     {
00346         coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
00347     }
00348     --order_;
00349 }
00350 
00351 template <class T>
00352 void
00353 PolynomialView<T>::backwardDeflate(T v)
00354 {
00355     v = -1.0 / v;
00356     coeffs_[0] *= v;
00357     for(unsigned int i = 1; i < order_; ++i)
00358     {
00359         coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
00360     }
00361     --order_;
00362 }
00363 
00364 template <class T>
00365 void
00366 PolynomialView<T>::deflateConjugatePair(Complex const & v)
00367 {
00368     vigra_precondition(order_ > 1,
00369         "PolynomialView<T>::deflateConjugatePair(): cannot deflate 2 roots "
00370         "from 1st order polynomial.");
00371     Real a = 2.0*v.real();
00372     Real b = -sq(v.real()) - sq(v.imag());
00373     coeffs_[order_-1] += a * coeffs_[order_];
00374     for(int i = order_-2; i > 1; --i)
00375     {
00376         coeffs_[i] += a * coeffs_[i+1] + b*coeffs_[i+2];
00377     }
00378     coeffs_ += 2;
00379     order_ -= 2;
00380 }
00381     
00382 template <class T>
00383 void 
00384 PolynomialView<T>::minimizeOrder(double epsilon)
00385 {
00386     while(std::abs(coeffs_[order_]) <= epsilon && order_ > 0)
00387             --order_;
00388 }
00389 
00390 template <class T>
00391 void 
00392 PolynomialView<T>::normalize()
00393 {
00394     for(unsigned int i = 0; i<order_; ++i)
00395         coeffs_[i] /= coeffs_[order_];
00396     coeffs_[order_] = T(1.0);
00397 }
00398 
00399 template <class T>
00400 void 
00401 PolynomialView<T>::balance()
00402 {
00403     Real p0 = abs(coeffs_[0]), po = abs(coeffs_[order_]);
00404     Real norm = (p0 > 0.0)
00405                     ? VIGRA_CSTD::sqrt(p0*po) 
00406                     : po;
00407     for(unsigned int i = 0; i<=order_; ++i)
00408         coeffs_[i] /= norm;
00409 }
00410 
00411 /*****************************************************************/
00412 /*                                                               */
00413 /*                           Polynomial                          */
00414 /*                                                               */
00415 /*****************************************************************/
00416 
00417 /** Polynomial with internally managed array.
00418 
00419     Most interesting functionality is inherited from \ref vigra::PolynomialView.
00420 
00421     \see vigra::PolynomialView, vigra::StaticPolynomial, polynomialRoots()
00422 
00423     <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br>
00424     Namespace: vigra
00425     
00426     \ingroup Polynomials
00427 */
00428 template <class T>
00429 class Polynomial
00430 : public PolynomialView<T>
00431 {
00432     typedef PolynomialView<T> BaseType;
00433   public:
00434     typedef typename BaseType::Real    Real;
00435     typedef typename BaseType::Complex Complex;
00436     typedef Polynomial<Real>           RealPolynomial;
00437     typedef Polynomial<Complex>        ComplexPolynomial;
00438 
00439     typedef T         value_type;
00440     typedef T *       iterator;
00441     typedef T const * const_iterator;    
00442     
00443         /** Construct polynomial with given <tt>order</tt> and all coefficients
00444             set to zero (they can be set later using <tt>operator[]</tt>
00445             or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 
00446             the precision of subsequent algorithms (especially root finding) 
00447             performed on the polynomial.
00448         */
00449     Polynomial(unsigned int order = 0, double epsilon = 1.0e-14)
00450     : BaseType(epsilon),
00451       polynomial_(order + 1, T())
00452     {
00453         setCoeffs(&polynomial_[0], order);
00454     }
00455     
00456         /** Copy constructor
00457         */
00458     Polynomial(Polynomial const & p)
00459     : BaseType(p.epsilon()),
00460       polynomial_(p.begin(), p.end())
00461     {
00462         setCoeffs(&polynomial_[0], p.order());
00463     }
00464 
00465         /** Construct polynomial by copying the given coefficient sequence.
00466         */
00467     template <class ITER>
00468     Polynomial(ITER i, unsigned int order)
00469     : BaseType(),
00470       polynomial_(i, i + order + 1)
00471     {
00472         setCoeffs(&polynomial_[0], order);
00473     }
00474     
00475         /** Construct polynomial by copying the given coefficient sequence.
00476             Set <tt>epsilon</tt> (default: 1.0e-14) as 
00477             the precision of subsequent algorithms (especially root finding) 
00478             performed on the polynomial.
00479         */
00480     template <class ITER>
00481     Polynomial(ITER i, unsigned int order, double epsilon)
00482     : BaseType(epsilon),
00483       polynomial_(i, i + order + 1)
00484     {
00485         setCoeffs(&polynomial_[0], order);
00486     }
00487     
00488         /** Assigment
00489         */
00490     Polynomial & operator=(Polynomial const & p)
00491     {
00492         if(this == &p)
00493             return *this;
00494         ArrayVector<T> tmp(p.begin(), p.end());
00495         polynomial_.swap(tmp);
00496         setCoeffs(&polynomial_[0], p.order());
00497         this->epsilon_ = p.epsilon_;
00498         return *this;
00499     }
00500     
00501         /** Construct new polynomial representing the derivative of this
00502             polynomial.
00503         */
00504     Polynomial<T> getDerivative(unsigned int n = 1) const
00505     {
00506         Polynomial<T> res(*this);
00507         res.differentiate(n);
00508         return res;
00509     }
00510 
00511         /** Construct new polynomial representing this polynomial after
00512             deflation at the real root <tt>r</tt>.
00513         */
00514     Polynomial<T> 
00515     getDeflated(Real r) const
00516     {
00517         Polynomial<T> res(*this);
00518         res.deflate(r);
00519         return res;
00520     }
00521 
00522         /** Construct new polynomial representing this polynomial after
00523             deflation at the complex root <tt>r</tt>. The resulting
00524             polynomial will have complex coefficients, even if this
00525             polynomial had real ones.
00526         */
00527     Polynomial<Complex> 
00528     getDeflated(Complex const & r) const
00529     {
00530         Polynomial<Complex> res(this->begin(), this->order(), this->epsilon());
00531         res.deflate(r);
00532         return res;
00533     }
00534 
00535   protected:
00536     ArrayVector<T> polynomial_;
00537 };
00538 
00539 /*****************************************************************/
00540 /*                                                               */
00541 /*                        StaticPolynomial                       */
00542 /*                                                               */
00543 /*****************************************************************/
00544 
00545 /** Polynomial with internally managed array of static length.
00546 
00547     Most interesting functionality is inherited from \ref vigra::PolynomialView.
00548     This class differs from \ref vigra::Polynomial in that it allocates
00549     its memory statically which is much faster. Therefore, <tt>StaticPolynomial</tt>
00550     can only represent polynomials up to the given <tt>MAXORDER</tt>.
00551 
00552     \see vigra::PolynomialView, vigra::Polynomial, polynomialRoots()
00553 
00554     <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br>
00555     Namespace: vigra
00556     
00557     \ingroup Polynomials
00558 */
00559 template <unsigned int MAXORDER, class T>
00560 class StaticPolynomial
00561 : public PolynomialView<T>
00562 {
00563     typedef PolynomialView<T> BaseType;
00564     
00565   public:
00566     typedef typename BaseType::Real    Real;
00567     typedef typename BaseType::Complex Complex;
00568     typedef StaticPolynomial<MAXORDER, Real> RealPolynomial;
00569     typedef StaticPolynomial<MAXORDER, Complex> ComplexPolynomial;
00570 
00571     typedef T         value_type;
00572     typedef T *       iterator;
00573     typedef T const * const_iterator;
00574     
00575     
00576         /** Construct polynomial with given <tt>order <= MAXORDER</tt> and all 
00577             coefficients set to zero (they can be set later using <tt>operator[]</tt>
00578             or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 
00579             the precision of subsequent algorithms (especially root finding) 
00580             performed on the polynomial.
00581         */
00582     StaticPolynomial(unsigned int order = 0, double epsilon = 1.0e-14)
00583     : BaseType(epsilon)
00584     {
00585         vigra_precondition(order <= MAXORDER,
00586             "StaticPolynomial(): order exceeds MAXORDER.");
00587         std::fill_n(polynomial_, order+1, T());
00588         setCoeffs(polynomial_, order);
00589     }
00590     
00591         /** Copy constructor
00592         */
00593     StaticPolynomial(StaticPolynomial const & p)
00594     : BaseType(p.epsilon())
00595     {
00596         std::copy(p.begin(), p.end(), polynomial_);
00597         setCoeffs(polynomial_, p.order());
00598     }
00599 
00600         /** Construct polynomial by copying the given coefficient sequence.
00601             <tt>order <= MAXORDER</tt> is required.
00602         */
00603     template <class ITER>
00604     StaticPolynomial(ITER i, unsigned int order)
00605     : BaseType()
00606     {
00607         vigra_precondition(order <= MAXORDER,
00608             "StaticPolynomial(): order exceeds MAXORDER.");
00609         std::copy(i, i + order + 1, polynomial_);
00610         setCoeffs(polynomial_, order);
00611     }
00612     
00613         /** Construct polynomial by copying the given coefficient sequence.
00614             <tt>order <= MAXORDER</tt> is required. Set <tt>epsilon</tt> (default: 1.0e-14) as 
00615             the precision of subsequent algorithms (especially root finding) 
00616             performed on the polynomial.
00617         */
00618     template <class ITER>
00619     StaticPolynomial(ITER i, unsigned int order, double epsilon)
00620     : BaseType(epsilon)
00621     {
00622         vigra_precondition(order <= MAXORDER,
00623             "StaticPolynomial(): order exceeds MAXORDER.");
00624         std::copy(i, i + order + 1, polynomial_);
00625         setCoeffs(polynomial_, order);
00626     }
00627     
00628         /** Assigment.
00629         */
00630     StaticPolynomial & operator=(StaticPolynomial const & p)
00631     {
00632         if(this == &p)
00633             return *this;
00634         std::copy(p.begin(), p.end(), polynomial_);
00635         setCoeffs(polynomial_, p.order());
00636         this->epsilon_ = p.epsilon_;
00637         return *this;
00638     }
00639     
00640         /** Construct new polynomial representing the derivative of this
00641             polynomial.
00642         */
00643     StaticPolynomial getDerivative(unsigned int n = 1) const
00644     {
00645         StaticPolynomial res(*this);
00646         res.differentiate(n);
00647         return res;
00648     }
00649 
00650         /** Construct new polynomial representing this polynomial after
00651             deflation at the real root <tt>r</tt>.
00652         */
00653     StaticPolynomial 
00654     getDeflated(Real r) const
00655     {
00656         StaticPolynomial res(*this);
00657         res.deflate(r);
00658         return res;
00659     }
00660 
00661         /** Construct new polynomial representing this polynomial after
00662             deflation at the complex root <tt>r</tt>. The resulting
00663             polynomial will have complex coefficients, even if this
00664             polynomial had real ones.
00665         */
00666     StaticPolynomial<MAXORDER, Complex> 
00667     getDeflated(Complex const & r) const
00668     {
00669         StaticPolynomial<MAXORDER, Complex>  res(this->begin(), this->order(), this->epsilon());
00670         res.deflate(r);
00671         return res;
00672     }
00673     
00674     void setOrder(unsigned int order)
00675     {
00676         vigra_precondition(order <= MAXORDER,
00677             "taticPolynomial::setOrder(): order exceeds MAXORDER.");
00678         this->order_ = order;
00679     }
00680 
00681   protected:
00682     T polynomial_[MAXORDER+1];
00683 };
00684 
00685 /************************************************************/
00686 
00687 namespace detail {
00688 
00689 // replacement for complex division (some compilers have numerically
00690 // less stable implementations); code from python complexobject.c
00691 template <class T>
00692 std::complex<T> complexDiv(std::complex<T> const & a, std::complex<T> const & b)
00693 {
00694      const double abs_breal = b.real() < 0 ? -b.real() : b.real();
00695      const double abs_bimag = b.imag() < 0 ? -b.imag() : b.imag();
00696 
00697      if (abs_breal >= abs_bimag) 
00698      {
00699         /* divide tops and bottom by b.real() */
00700         if (abs_breal == 0.0) 
00701         {
00702             return std::complex<T>(a.real() / abs_breal, a.imag() / abs_breal);
00703         }
00704         else 
00705         {
00706             const double ratio = b.imag() / b.real();
00707             const double denom = b.real() + b.imag() * ratio;
00708             return std::complex<T>((a.real() + a.imag() * ratio) / denom,
00709                                    (a.imag() - a.real() * ratio) / denom);
00710         }
00711     }
00712     else 
00713     {
00714         /* divide tops and bottom by b.imag() */
00715         const double ratio = b.real() / b.imag();
00716         const double denom = b.real() * ratio + b.imag();
00717         return std::complex<T>((a.real() * ratio + a.imag()) / denom,
00718                                (a.imag() * ratio - a.real()) / denom);
00719     }
00720 }
00721 
00722 template <class T>
00723 std::complex<T> deleteBelowEpsilon(std::complex<T> const & x, double eps)
00724 {
00725     return std::abs(x.imag()) <= 2.0*eps*std::abs(x.real()) 
00726                 ? std::complex<T>(x.real())
00727                 : std::abs(x.real()) <= 2.0*eps*std::abs(x.imag()) 
00728                     ? std::complex<T>(NumericTraits<T>::zero(), x.imag())
00729                     :  x;
00730 }
00731 
00732 template <class POLYNOMIAL>
00733 typename POLYNOMIAL::value_type
00734 laguerreStartingGuess(POLYNOMIAL const & p)
00735 {
00736     double N = p.order();
00737     typename POLYNOMIAL::value_type centroid = -p[p.order()-1] / N / p[p.order()];
00738     double dist = VIGRA_CSTD::pow(std::abs(p(centroid) / p[p.order()]), 1.0 / N);
00739     return centroid + dist;
00740 }
00741 
00742 template <class POLYNOMIAL, class Complex>
00743 int laguerre1Root(POLYNOMIAL const & p, Complex & x, unsigned int multiplicity)
00744 {
00745     typedef typename NumericTraits<Complex>::ValueType Real;
00746     
00747     static double frac[] = {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
00748     int maxiter = 80, 
00749         count;
00750     double N = p.order();
00751     double eps  = p.epsilon(),
00752            eps2 = VIGRA_CSTD::sqrt(eps);
00753         
00754     if(multiplicity == 0)
00755         x = laguerreStartingGuess(p);
00756         
00757     bool mayTryDerivative = true;  // try derivative for multiple roots
00758     
00759     for(count = 0; count < maxiter; ++count)
00760     {
00761         // Horner's algorithm to calculate values of polynomial and its
00762         // first two derivatives and estimate error for current x
00763         Complex p0(p[p.order()]);
00764         Complex p1(0.0);
00765         Complex p2(0.0);
00766         Real ax    = std::abs(x);
00767         Real err = std::abs(p0);
00768         for(int i = p.order()-1; i >= 0; --i)
00769         {
00770             p2  = p2  * x  + p1;
00771             p1  = p1  * x  + p0;
00772             p0  = p0  * x  + p[i];
00773             err = err * ax + std::abs(p0);
00774         }
00775         p2 *= 2.0;
00776         err *= eps;
00777         Real ap0 = std::abs(p0);
00778         if(ap0 <= err)
00779         {
00780             break;  // converged
00781         }
00782         Complex g = complexDiv(p1, p0);
00783         Complex g2 = g * g;
00784         Complex h = g2 - complexDiv(p2, p0);
00785         // estimate root multiplicity according to Tien Chen
00786         if(g2 != 0.0)
00787         {
00788             multiplicity = (unsigned int)VIGRA_CSTD::floor(N / 
00789                                 (std::abs(N * complexDiv(h, g2) - 1.0) + 1.0) + 0.5);
00790             if(multiplicity < 1)
00791                 multiplicity = 1;
00792         }
00793         // improve accuracy of multiple roots on the derivative, as suggested by C. Bond
00794         // (do this only if we are already near the root, otherwise we may converge to 
00795         //  a different root of the derivative polynomial)
00796         if(mayTryDerivative && multiplicity > 1 && ap0 < eps2)
00797         {
00798             Complex x1 = x;
00799             int derivativeMultiplicity = laguerre1Root(p.getDerivative(), x1, multiplicity-1);
00800             if(derivativeMultiplicity && std::abs(p(x1)) < std::abs(p(x)))
00801             {
00802                 // successful search on derivative
00803                 x = x1;
00804                 return derivativeMultiplicity + 1;
00805             }
00806             else
00807             {
00808                 // unsuccessful search on derivative => don't do it again
00809                 mayTryDerivative = false;
00810             }
00811         }
00812         Complex sq = VIGRA_CSTD::sqrt((N - 1.0) * (N * h - g2));
00813         Complex gp = g + sq;
00814         Complex gm = g - sq;
00815         if(std::abs(gp) < std::abs(gm))
00816             gp = gm;
00817         Complex dx;
00818         if(gp != 0.0)
00819         {
00820             dx = complexDiv(Complex(N) , gp);
00821         }
00822         else
00823         {
00824             // re-initialisation trick due to Numerical Recipes
00825             dx = (1.0 + ax) * Complex(VIGRA_CSTD::cos(double(count)), VIGRA_CSTD::sin(double(count)));
00826         }
00827         Complex x1 = x - dx;
00828 
00829         if(x1 - x == 0.0)
00830         {
00831             break;  // converged
00832         }
00833         if((count + 1) % 10)
00834             x = x1;
00835         else
00836             // cycle breaking trick according to Numerical Recipes
00837             x = x - frac[(count+1)/10] * dx;
00838     }
00839     return count < maxiter ? 
00840         multiplicity : 
00841         0;
00842 }
00843 
00844 template <class Real>
00845 struct PolynomialRootCompare
00846 {
00847     Real epsilon;
00848 
00849     PolynomialRootCompare(Real eps)
00850     : epsilon(eps)
00851     {}
00852     
00853     template <class T>
00854     bool operator()(T const & l, T const & r)
00855     {
00856         return closeAtTolerance(l.real(), r.real(), epsilon)
00857                      ? l.imag() < r.imag()
00858                      : l.real() < r.real();
00859     }
00860 };
00861 
00862 } // namespace detail 
00863 
00864 /** \addtogroup Polynomials Polynomials and root determination
00865 
00866     Classes to represent polynomials and functions to find polynomial roots.
00867 */
00868 //@{
00869 
00870 /*****************************************************************/
00871 /*                                                               */
00872 /*                         polynomialRoots                       */
00873 /*                                                               */
00874 /*****************************************************************/
00875 
00876 /** Determine the roots of the polynomial <tt>poriginal</tt>.
00877 
00878     The roots are appended to the vector <tt>roots</tt>, with optional root
00879     polishing as specified by <tt>polishRoots</tt> (default: do polishing). The function uses an 
00880     improved version of Laguerre's algorithm. The improvements are as follows:
00881     
00882     <ul>
00883     <li>It uses a clever initial guess for the iteration, according to a proposal by Tien Chen</li>
00884     <li>It estimates each root's multiplicity, again according to Tien Chen, and reduces multiplicity
00885         by switching to the polynomial's derivative (which has the same root, with multiplicity
00886         reduced by one), as proposed by C. Bond.</li>
00887     </ul>
00888     
00889     The algorithm has been successfully used for polynomials up to order 80.
00890     The function stops and returns <tt>false</tt> if an iteration fails to converge within 
00891     80 steps. The type <tt>POLYNOMIAL</tt> must be compatible to 
00892     \ref vigra::PolynomialView, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
00893     with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>.
00894 
00895     <b> Declaration:</b>
00896 
00897     pass arguments explicitly:
00898     \code
00899     namespace vigra {
00900         template <class POLYNOMIAL, class VECTOR>
00901         bool 
00902         polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots = true);
00903     }
00904     \endcode
00905 
00906 
00907     <b> Usage:</b>
00908 
00909         <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br>
00910         Namespace: vigra
00911 
00912     \code
00913     // encode the polynomial  x^4 - 1
00914     Polynomial<double> poly(4);
00915     poly[0] = -1.0;
00916     poly[4] =  1.0;
00917 
00918     ArrayVector<std::complex<double> > roots;
00919     polynomialRoots(poly, roots);
00920     \endcode
00921 
00922     \see polynomialRootsEigenvalueMethod()
00923 */
00924 template <class POLYNOMIAL, class VECTOR>
00925 bool polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots)
00926 {
00927     typedef typename POLYNOMIAL::value_type T;
00928     typedef typename POLYNOMIAL::Real    Real;
00929     typedef typename POLYNOMIAL::Complex Complex;
00930     typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial;
00931     
00932     double eps  = poriginal.epsilon();
00933 
00934     WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
00935     p.minimizeOrder();
00936     if(p.order() == 0)
00937         return true;
00938 
00939     Complex x = detail::laguerreStartingGuess(p);
00940     
00941     unsigned int multiplicity = 1;
00942     bool triedConjugate = false;
00943     
00944     // handle the high order cases
00945     while(p.order() > 2)
00946     {
00947         p.balance();
00948         
00949         // find root estimate using Laguerre's method on deflated polynomial p;
00950         // zero return indicates failure to converge
00951         multiplicity = detail::laguerre1Root(p, x, multiplicity);
00952     
00953         if(multiplicity == 0)
00954             return false;
00955         // polish root on original polynomial poriginal;
00956         // zero return indicates failure to converge
00957         if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
00958             return false;
00959         x = detail::deleteBelowEpsilon(x, eps);
00960         roots.push_back(x);
00961         p.deflate(x);
00962         // determine the next starting guess
00963         if(multiplicity > 1)
00964         {
00965             // probably multiple root => keep current root as starting guess
00966             --multiplicity;
00967             triedConjugate = false;
00968         }
00969         else
00970         {
00971             // need a new starting guess
00972             if(x.imag() != 0.0 && !triedConjugate)
00973             {
00974                 // if the root is complex and we don't already have 
00975                 // the conjugate root => try the conjugate as starting guess
00976                 triedConjugate = true;
00977                 x = conj(x);
00978             }
00979             else
00980             {
00981                 // otherwise generate new starting guess
00982                 triedConjugate = false;
00983                 x = detail::laguerreStartingGuess(p);
00984             }
00985         }
00986     }
00987     
00988     // handle the low order cases
00989     if(p.order() == 2)
00990     {
00991         Complex a = p[2];
00992         Complex b = p[1];
00993         Complex c = p[0];
00994         Complex b2 = std::sqrt(b*b - 4.0*a*c);
00995         Complex q;
00996         if((conj(b)*b2).real() >= 0.0)
00997             q = -0.5 * (b + b2);
00998         else
00999             q = -0.5 * (b - b2);
01000         x = detail::complexDiv(q, a);
01001         if(polishRoots)
01002             detail::laguerre1Root(poriginal, x, 1);
01003         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01004         x = detail::complexDiv(c, q);
01005         if(polishRoots)
01006             detail::laguerre1Root(poriginal, x, 1);
01007         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01008     }
01009     else if(p.order() == 1)
01010     {
01011         x = detail::complexDiv(-p[0], p[1]);
01012         if(polishRoots)
01013             detail::laguerre1Root(poriginal, x, 1);
01014         roots.push_back(detail::deleteBelowEpsilon(x, eps));
01015     }
01016     std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
01017     return true;
01018 }
01019 
01020 template <class POLYNOMIAL, class VECTOR>
01021 inline bool 
01022 polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
01023 {
01024     return polynomialRoots(poriginal, roots, true);
01025 }
01026 
01027 /** Determine the real roots of the polynomial <tt>p</tt>.
01028 
01029     This function simply calls \ref polynomialRoots() and than throws away all complex roots.
01030     Accordingly, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
01031     with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>.
01032 
01033     <b> Declaration:</b>
01034 
01035     pass arguments explicitly:
01036     \code
01037     namespace vigra {
01038         template <class POLYNOMIAL, class VECTOR>
01039         bool 
01040         polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots = true);
01041     }
01042     \endcode
01043 
01044 
01045     <b> Usage:</b>
01046 
01047         <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br>
01048         Namespace: vigra
01049 
01050     \code
01051     // encode the polynomial  x^4 - 1
01052     Polynomial<double> poly(4);
01053     poly[0] = -1.0;
01054     poly[4] =  1.0;
01055 
01056     ArrayVector<double> roots;
01057     polynomialRealRoots(poly, roots);
01058     \endcode
01059 
01060     \see polynomialRealRootsEigenvalueMethod()
01061 */
01062 template <class POLYNOMIAL, class VECTOR>
01063 bool polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots)
01064 {
01065     typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
01066     ArrayVector<Complex> croots;
01067     if(!polynomialRoots(p, croots, polishRoots))
01068         return false;
01069     for(unsigned int i = 0; i < croots.size(); ++i)
01070         if(croots[i].imag() == 0.0)
01071             roots.push_back(croots[i].real());
01072     return true;
01073 }
01074 
01075 template <class POLYNOMIAL, class VECTOR>
01076 inline bool 
01077 polynomialRealRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
01078 {
01079     return polynomialRealRoots(poriginal, roots, true);
01080 }
01081 
01082 //@}
01083 
01084 } // namespace vigra
01085 
01086 #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.3 (18 Aug 2005)