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