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