[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><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 ontermediate 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 snallest 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> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><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> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><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> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><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> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><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) |
html generated using doxygen and Python
|