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