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

vigra/polynomial.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)