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

vigra/polynomial.hxx

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

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

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)