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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)