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