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

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