[ 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.3, Aug 18 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 void balance(); 00188 00189 /** Get iterator for the coefficient sequence. 00190 */ 00191 iterator begin() 00192 { return coeffs_; } 00193 00194 /** Get end iterator for the coefficient sequence. 00195 */ 00196 iterator end() 00197 { return begin() + size(); } 00198 00199 /** Get const_iterator for the coefficient sequence. 00200 */ 00201 const_iterator begin() const 00202 { return coeffs_; } 00203 00204 /** Get end const_iterator for the coefficient sequence. 00205 */ 00206 const_iterator end() const 00207 { return begin() + size(); } 00208 00209 /** Get length of the coefficient sequence (<tt>order() + 1</tt>). 00210 */ 00211 unsigned int size() const 00212 { return order_ + 1; } 00213 00214 /** Get order of the polynomial. 00215 */ 00216 unsigned int order() const 00217 { return order_; } 00218 00219 /** Get requested precision for polynomial algorithms 00220 (especially root finding). 00221 */ 00222 double epsilon() const 00223 { return epsilon_; } 00224 00225 /** Set requested precision for polynomial algorithms 00226 (especially root finding). 00227 */ 00228 void setEpsilon(double eps) 00229 { epsilon_ = eps; } 00230 00231 protected: 00232 PolynomialView(double epsilon = 1e-14) 00233 : coeffs_(0), 00234 order_(0), 00235 epsilon_(epsilon) 00236 {} 00237 00238 void setCoeffs(T * coeffs, unsigned int order) 00239 { 00240 coeffs_ = coeffs; 00241 order_ = order; 00242 } 00243 00244 T * coeffs_; 00245 unsigned int order_; 00246 double epsilon_; 00247 }; 00248 00249 template <class T> 00250 template <class U> 00251 typename PromoteTraits<T, U>::Promote 00252 PolynomialView<T>::operator()(U v) const 00253 { 00254 typename PromoteTraits<T, U>::Promote p(coeffs_[order_]); 00255 for(int i = order_ - 1; i >= 0; --i) 00256 { 00257 p = v * p + coeffs_[i]; 00258 } 00259 return p; 00260 } 00261 00262 /* 00263 template <class T> 00264 typename PolynomialView<T>::Complex 00265 PolynomialView<T>::operator()(Complex const & v) const 00266 { 00267 Complex p(coeffs_[order_]); 00268 for(int i = order_ - 1; i >= 0; --i) 00269 { 00270 p = v * p + coeffs_[i]; 00271 } 00272 return p; 00273 } 00274 */ 00275 00276 template <class T> 00277 void 00278 PolynomialView<T>::differentiate(unsigned int n) 00279 { 00280 if(n == 0) 00281 return; 00282 if(order_ == 0) 00283 { 00284 coeffs_[0] = 0.0; 00285 return; 00286 } 00287 for(unsigned int i = 1; i <= order_; ++i) 00288 { 00289 coeffs_[i-1] = double(i)*coeffs_[i]; 00290 } 00291 --order_; 00292 if(n > 1) 00293 differentiate(n-1); 00294 } 00295 00296 template <class T> 00297 void 00298 PolynomialView<T>::deflate(T const & v, unsigned int multiplicity) 00299 { 00300 vigra_precondition(order_ > 0, 00301 "PolynomialView<T>::deflate(): cannot deflate 0th order polynomial."); 00302 if(v == 0.0) 00303 { 00304 ++coeffs_; 00305 --order_; 00306 } 00307 else 00308 { 00309 // we use combined forward/backward deflation because 00310 // our initial guess seems to favour convergence to 00311 // a root with magnitude near the median among all roots 00312 forwardBackwardDeflate(v); 00313 } 00314 if(multiplicity > 1) 00315 deflate(v, multiplicity-1); 00316 } 00317 00318 template <class T> 00319 void 00320 PolynomialView<T>::forwardDeflate(T const & v) 00321 { 00322 for(int i = order_-1; i > 0; --i) 00323 { 00324 coeffs_[i] += v * coeffs_[i+1]; 00325 } 00326 ++coeffs_; 00327 --order_; 00328 } 00329 00330 template <class T> 00331 void 00332 PolynomialView<T>::forwardBackwardDeflate(T v) 00333 { 00334 unsigned int order2 = order_ / 2; 00335 T tmp = coeffs_[order_]; 00336 for(unsigned int i = order_-1; i >= order2; --i) 00337 { 00338 T tmp1 = coeffs_[i]; 00339 coeffs_[i] = tmp; 00340 tmp = tmp1 + v * tmp; 00341 } 00342 v = -1.0 / v; 00343 coeffs_[0] *= v; 00344 for(unsigned int i = 1; i < order2; ++i) 00345 { 00346 coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]); 00347 } 00348 --order_; 00349 } 00350 00351 template <class T> 00352 void 00353 PolynomialView<T>::backwardDeflate(T v) 00354 { 00355 v = -1.0 / v; 00356 coeffs_[0] *= v; 00357 for(unsigned int i = 1; i < order_; ++i) 00358 { 00359 coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]); 00360 } 00361 --order_; 00362 } 00363 00364 template <class T> 00365 void 00366 PolynomialView<T>::deflateConjugatePair(Complex const & v) 00367 { 00368 vigra_precondition(order_ > 1, 00369 "PolynomialView<T>::deflateConjugatePair(): cannot deflate 2 roots " 00370 "from 1st order polynomial."); 00371 Real a = 2.0*v.real(); 00372 Real b = -sq(v.real()) - sq(v.imag()); 00373 coeffs_[order_-1] += a * coeffs_[order_]; 00374 for(int i = order_-2; i > 1; --i) 00375 { 00376 coeffs_[i] += a * coeffs_[i+1] + b*coeffs_[i+2]; 00377 } 00378 coeffs_ += 2; 00379 order_ -= 2; 00380 } 00381 00382 template <class T> 00383 void 00384 PolynomialView<T>::minimizeOrder(double epsilon) 00385 { 00386 while(std::abs(coeffs_[order_]) <= epsilon && order_ > 0) 00387 --order_; 00388 } 00389 00390 template <class T> 00391 void 00392 PolynomialView<T>::normalize() 00393 { 00394 for(unsigned int i = 0; i<order_; ++i) 00395 coeffs_[i] /= coeffs_[order_]; 00396 coeffs_[order_] = T(1.0); 00397 } 00398 00399 template <class T> 00400 void 00401 PolynomialView<T>::balance() 00402 { 00403 Real p0 = abs(coeffs_[0]), po = abs(coeffs_[order_]); 00404 Real norm = (p0 > 0.0) 00405 ? VIGRA_CSTD::sqrt(p0*po) 00406 : po; 00407 for(unsigned int i = 0; i<=order_; ++i) 00408 coeffs_[i] /= norm; 00409 } 00410 00411 /*****************************************************************/ 00412 /* */ 00413 /* Polynomial */ 00414 /* */ 00415 /*****************************************************************/ 00416 00417 /** Polynomial with internally managed array. 00418 00419 Most interesting functionality is inherited from \ref vigra::PolynomialView. 00420 00421 \see vigra::PolynomialView, vigra::StaticPolynomial, polynomialRoots() 00422 00423 <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br> 00424 Namespace: vigra 00425 00426 \ingroup Polynomials 00427 */ 00428 template <class T> 00429 class Polynomial 00430 : public PolynomialView<T> 00431 { 00432 typedef PolynomialView<T> BaseType; 00433 public: 00434 typedef typename BaseType::Real Real; 00435 typedef typename BaseType::Complex Complex; 00436 typedef Polynomial<Real> RealPolynomial; 00437 typedef Polynomial<Complex> ComplexPolynomial; 00438 00439 typedef T value_type; 00440 typedef T * iterator; 00441 typedef T const * const_iterator; 00442 00443 /** Construct polynomial with given <tt>order</tt> and all coefficients 00444 set to zero (they can be set later using <tt>operator[]</tt> 00445 or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 00446 the precision of subsequent algorithms (especially root finding) 00447 performed on the polynomial. 00448 */ 00449 Polynomial(unsigned int order = 0, double epsilon = 1.0e-14) 00450 : BaseType(epsilon), 00451 polynomial_(order + 1, T()) 00452 { 00453 setCoeffs(&polynomial_[0], order); 00454 } 00455 00456 /** Copy constructor 00457 */ 00458 Polynomial(Polynomial const & p) 00459 : BaseType(p.epsilon()), 00460 polynomial_(p.begin(), p.end()) 00461 { 00462 setCoeffs(&polynomial_[0], p.order()); 00463 } 00464 00465 /** Construct polynomial by copying the given coefficient sequence. 00466 */ 00467 template <class ITER> 00468 Polynomial(ITER i, unsigned int order) 00469 : BaseType(), 00470 polynomial_(i, i + order + 1) 00471 { 00472 setCoeffs(&polynomial_[0], order); 00473 } 00474 00475 /** Construct polynomial by copying the given coefficient sequence. 00476 Set <tt>epsilon</tt> (default: 1.0e-14) as 00477 the precision of subsequent algorithms (especially root finding) 00478 performed on the polynomial. 00479 */ 00480 template <class ITER> 00481 Polynomial(ITER i, unsigned int order, double epsilon) 00482 : BaseType(epsilon), 00483 polynomial_(i, i + order + 1) 00484 { 00485 setCoeffs(&polynomial_[0], order); 00486 } 00487 00488 /** Assigment 00489 */ 00490 Polynomial & operator=(Polynomial const & p) 00491 { 00492 if(this == &p) 00493 return *this; 00494 ArrayVector<T> tmp(p.begin(), p.end()); 00495 polynomial_.swap(tmp); 00496 setCoeffs(&polynomial_[0], p.order()); 00497 this->epsilon_ = p.epsilon_; 00498 return *this; 00499 } 00500 00501 /** Construct new polynomial representing the derivative of this 00502 polynomial. 00503 */ 00504 Polynomial<T> getDerivative(unsigned int n = 1) const 00505 { 00506 Polynomial<T> res(*this); 00507 res.differentiate(n); 00508 return res; 00509 } 00510 00511 /** Construct new polynomial representing this polynomial after 00512 deflation at the real root <tt>r</tt>. 00513 */ 00514 Polynomial<T> 00515 getDeflated(Real r) const 00516 { 00517 Polynomial<T> res(*this); 00518 res.deflate(r); 00519 return res; 00520 } 00521 00522 /** Construct new polynomial representing this polynomial after 00523 deflation at the complex root <tt>r</tt>. The resulting 00524 polynomial will have complex coefficients, even if this 00525 polynomial had real ones. 00526 */ 00527 Polynomial<Complex> 00528 getDeflated(Complex const & r) const 00529 { 00530 Polynomial<Complex> res(this->begin(), this->order(), this->epsilon()); 00531 res.deflate(r); 00532 return res; 00533 } 00534 00535 protected: 00536 ArrayVector<T> polynomial_; 00537 }; 00538 00539 /*****************************************************************/ 00540 /* */ 00541 /* StaticPolynomial */ 00542 /* */ 00543 /*****************************************************************/ 00544 00545 /** Polynomial with internally managed array of static length. 00546 00547 Most interesting functionality is inherited from \ref vigra::PolynomialView. 00548 This class differs from \ref vigra::Polynomial in that it allocates 00549 its memory statically which is much faster. Therefore, <tt>StaticPolynomial</tt> 00550 can only represent polynomials up to the given <tt>MAXORDER</tt>. 00551 00552 \see vigra::PolynomialView, vigra::Polynomial, polynomialRoots() 00553 00554 <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br> 00555 Namespace: vigra 00556 00557 \ingroup Polynomials 00558 */ 00559 template <unsigned int MAXORDER, class T> 00560 class StaticPolynomial 00561 : public PolynomialView<T> 00562 { 00563 typedef PolynomialView<T> BaseType; 00564 00565 public: 00566 typedef typename BaseType::Real Real; 00567 typedef typename BaseType::Complex Complex; 00568 typedef StaticPolynomial<MAXORDER, Real> RealPolynomial; 00569 typedef StaticPolynomial<MAXORDER, Complex> ComplexPolynomial; 00570 00571 typedef T value_type; 00572 typedef T * iterator; 00573 typedef T const * const_iterator; 00574 00575 00576 /** Construct polynomial with given <tt>order <= MAXORDER</tt> and all 00577 coefficients set to zero (they can be set later using <tt>operator[]</tt> 00578 or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines 00579 the precision of subsequent algorithms (especially root finding) 00580 performed on the polynomial. 00581 */ 00582 StaticPolynomial(unsigned int order = 0, double epsilon = 1.0e-14) 00583 : BaseType(epsilon) 00584 { 00585 vigra_precondition(order <= MAXORDER, 00586 "StaticPolynomial(): order exceeds MAXORDER."); 00587 std::fill_n(polynomial_, order+1, T()); 00588 setCoeffs(polynomial_, order); 00589 } 00590 00591 /** Copy constructor 00592 */ 00593 StaticPolynomial(StaticPolynomial const & p) 00594 : BaseType(p.epsilon()) 00595 { 00596 std::copy(p.begin(), p.end(), polynomial_); 00597 setCoeffs(polynomial_, p.order()); 00598 } 00599 00600 /** Construct polynomial by copying the given coefficient sequence. 00601 <tt>order <= MAXORDER</tt> is required. 00602 */ 00603 template <class ITER> 00604 StaticPolynomial(ITER i, unsigned int order) 00605 : BaseType() 00606 { 00607 vigra_precondition(order <= MAXORDER, 00608 "StaticPolynomial(): order exceeds MAXORDER."); 00609 std::copy(i, i + order + 1, polynomial_); 00610 setCoeffs(polynomial_, order); 00611 } 00612 00613 /** Construct polynomial by copying the given coefficient sequence. 00614 <tt>order <= MAXORDER</tt> is required. Set <tt>epsilon</tt> (default: 1.0e-14) as 00615 the precision of subsequent algorithms (especially root finding) 00616 performed on the polynomial. 00617 */ 00618 template <class ITER> 00619 StaticPolynomial(ITER i, unsigned int order, double epsilon) 00620 : BaseType(epsilon) 00621 { 00622 vigra_precondition(order <= MAXORDER, 00623 "StaticPolynomial(): order exceeds MAXORDER."); 00624 std::copy(i, i + order + 1, polynomial_); 00625 setCoeffs(polynomial_, order); 00626 } 00627 00628 /** Assigment. 00629 */ 00630 StaticPolynomial & operator=(StaticPolynomial const & p) 00631 { 00632 if(this == &p) 00633 return *this; 00634 std::copy(p.begin(), p.end(), polynomial_); 00635 setCoeffs(polynomial_, p.order()); 00636 this->epsilon_ = p.epsilon_; 00637 return *this; 00638 } 00639 00640 /** Construct new polynomial representing the derivative of this 00641 polynomial. 00642 */ 00643 StaticPolynomial getDerivative(unsigned int n = 1) const 00644 { 00645 StaticPolynomial res(*this); 00646 res.differentiate(n); 00647 return res; 00648 } 00649 00650 /** Construct new polynomial representing this polynomial after 00651 deflation at the real root <tt>r</tt>. 00652 */ 00653 StaticPolynomial 00654 getDeflated(Real r) const 00655 { 00656 StaticPolynomial res(*this); 00657 res.deflate(r); 00658 return res; 00659 } 00660 00661 /** Construct new polynomial representing this polynomial after 00662 deflation at the complex root <tt>r</tt>. The resulting 00663 polynomial will have complex coefficients, even if this 00664 polynomial had real ones. 00665 */ 00666 StaticPolynomial<MAXORDER, Complex> 00667 getDeflated(Complex const & r) const 00668 { 00669 StaticPolynomial<MAXORDER, Complex> res(this->begin(), this->order(), this->epsilon()); 00670 res.deflate(r); 00671 return res; 00672 } 00673 00674 void setOrder(unsigned int order) 00675 { 00676 vigra_precondition(order <= MAXORDER, 00677 "taticPolynomial::setOrder(): order exceeds MAXORDER."); 00678 this->order_ = order; 00679 } 00680 00681 protected: 00682 T polynomial_[MAXORDER+1]; 00683 }; 00684 00685 /************************************************************/ 00686 00687 namespace detail { 00688 00689 // replacement for complex division (some compilers have numerically 00690 // less stable implementations); code from python complexobject.c 00691 template <class T> 00692 std::complex<T> complexDiv(std::complex<T> const & a, std::complex<T> const & b) 00693 { 00694 const double abs_breal = b.real() < 0 ? -b.real() : b.real(); 00695 const double abs_bimag = b.imag() < 0 ? -b.imag() : b.imag(); 00696 00697 if (abs_breal >= abs_bimag) 00698 { 00699 /* divide tops and bottom by b.real() */ 00700 if (abs_breal == 0.0) 00701 { 00702 return std::complex<T>(a.real() / abs_breal, a.imag() / abs_breal); 00703 } 00704 else 00705 { 00706 const double ratio = b.imag() / b.real(); 00707 const double denom = b.real() + b.imag() * ratio; 00708 return std::complex<T>((a.real() + a.imag() * ratio) / denom, 00709 (a.imag() - a.real() * ratio) / denom); 00710 } 00711 } 00712 else 00713 { 00714 /* divide tops and bottom by b.imag() */ 00715 const double ratio = b.real() / b.imag(); 00716 const double denom = b.real() * ratio + b.imag(); 00717 return std::complex<T>((a.real() * ratio + a.imag()) / denom, 00718 (a.imag() * ratio - a.real()) / denom); 00719 } 00720 } 00721 00722 template <class T> 00723 std::complex<T> deleteBelowEpsilon(std::complex<T> const & x, double eps) 00724 { 00725 return std::abs(x.imag()) <= 2.0*eps*std::abs(x.real()) 00726 ? std::complex<T>(x.real()) 00727 : std::abs(x.real()) <= 2.0*eps*std::abs(x.imag()) 00728 ? std::complex<T>(NumericTraits<T>::zero(), x.imag()) 00729 : x; 00730 } 00731 00732 template <class POLYNOMIAL> 00733 typename POLYNOMIAL::value_type 00734 laguerreStartingGuess(POLYNOMIAL const & p) 00735 { 00736 double N = p.order(); 00737 typename POLYNOMIAL::value_type centroid = -p[p.order()-1] / N / p[p.order()]; 00738 double dist = VIGRA_CSTD::pow(std::abs(p(centroid) / p[p.order()]), 1.0 / N); 00739 return centroid + dist; 00740 } 00741 00742 template <class POLYNOMIAL, class Complex> 00743 int laguerre1Root(POLYNOMIAL const & p, Complex & x, unsigned int multiplicity) 00744 { 00745 typedef typename NumericTraits<Complex>::ValueType Real; 00746 00747 static double frac[] = {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0}; 00748 int maxiter = 80, 00749 count; 00750 double N = p.order(); 00751 double eps = p.epsilon(), 00752 eps2 = VIGRA_CSTD::sqrt(eps); 00753 00754 if(multiplicity == 0) 00755 x = laguerreStartingGuess(p); 00756 00757 bool mayTryDerivative = true; // try derivative for multiple roots 00758 00759 for(count = 0; count < maxiter; ++count) 00760 { 00761 // Horner's algorithm to calculate values of polynomial and its 00762 // first two derivatives and estimate error for current x 00763 Complex p0(p[p.order()]); 00764 Complex p1(0.0); 00765 Complex p2(0.0); 00766 Real ax = std::abs(x); 00767 Real err = std::abs(p0); 00768 for(int i = p.order()-1; i >= 0; --i) 00769 { 00770 p2 = p2 * x + p1; 00771 p1 = p1 * x + p0; 00772 p0 = p0 * x + p[i]; 00773 err = err * ax + std::abs(p0); 00774 } 00775 p2 *= 2.0; 00776 err *= eps; 00777 Real ap0 = std::abs(p0); 00778 if(ap0 <= err) 00779 { 00780 break; // converged 00781 } 00782 Complex g = complexDiv(p1, p0); 00783 Complex g2 = g * g; 00784 Complex h = g2 - complexDiv(p2, p0); 00785 // estimate root multiplicity according to Tien Chen 00786 if(g2 != 0.0) 00787 { 00788 multiplicity = (unsigned int)VIGRA_CSTD::floor(N / 00789 (std::abs(N * complexDiv(h, g2) - 1.0) + 1.0) + 0.5); 00790 if(multiplicity < 1) 00791 multiplicity = 1; 00792 } 00793 // improve accuracy of multiple roots on the derivative, as suggested by C. Bond 00794 // (do this only if we are already near the root, otherwise we may converge to 00795 // a different root of the derivative polynomial) 00796 if(mayTryDerivative && multiplicity > 1 && ap0 < eps2) 00797 { 00798 Complex x1 = x; 00799 int derivativeMultiplicity = laguerre1Root(p.getDerivative(), x1, multiplicity-1); 00800 if(derivativeMultiplicity && std::abs(p(x1)) < std::abs(p(x))) 00801 { 00802 // successful search on derivative 00803 x = x1; 00804 return derivativeMultiplicity + 1; 00805 } 00806 else 00807 { 00808 // unsuccessful search on derivative => don't do it again 00809 mayTryDerivative = false; 00810 } 00811 } 00812 Complex sq = VIGRA_CSTD::sqrt((N - 1.0) * (N * h - g2)); 00813 Complex gp = g + sq; 00814 Complex gm = g - sq; 00815 if(std::abs(gp) < std::abs(gm)) 00816 gp = gm; 00817 Complex dx; 00818 if(gp != 0.0) 00819 { 00820 dx = complexDiv(Complex(N) , gp); 00821 } 00822 else 00823 { 00824 // re-initialisation trick due to Numerical Recipes 00825 dx = (1.0 + ax) * Complex(VIGRA_CSTD::cos(double(count)), VIGRA_CSTD::sin(double(count))); 00826 } 00827 Complex x1 = x - dx; 00828 00829 if(x1 - x == 0.0) 00830 { 00831 break; // converged 00832 } 00833 if((count + 1) % 10) 00834 x = x1; 00835 else 00836 // cycle breaking trick according to Numerical Recipes 00837 x = x - frac[(count+1)/10] * dx; 00838 } 00839 return count < maxiter ? 00840 multiplicity : 00841 0; 00842 } 00843 00844 template <class Real> 00845 struct PolynomialRootCompare 00846 { 00847 Real epsilon; 00848 00849 PolynomialRootCompare(Real eps) 00850 : epsilon(eps) 00851 {} 00852 00853 template <class T> 00854 bool operator()(T const & l, T const & r) 00855 { 00856 return closeAtTolerance(l.real(), r.real(), epsilon) 00857 ? l.imag() < r.imag() 00858 : l.real() < r.real(); 00859 } 00860 }; 00861 00862 } // namespace detail 00863 00864 /** \addtogroup Polynomials Polynomials and root determination 00865 00866 Classes to represent polynomials and functions to find polynomial roots. 00867 */ 00868 //@{ 00869 00870 /*****************************************************************/ 00871 /* */ 00872 /* polynomialRoots */ 00873 /* */ 00874 /*****************************************************************/ 00875 00876 /** Determine the roots of the polynomial <tt>poriginal</tt>. 00877 00878 The roots are appended to the vector <tt>roots</tt>, with optional root 00879 polishing as specified by <tt>polishRoots</tt> (default: do polishing). The function uses an 00880 improved version of Laguerre's algorithm. The improvements are as follows: 00881 00882 <ul> 00883 <li>It uses a clever initial guess for the iteration, according to a proposal by Tien Chen</li> 00884 <li>It estimates each root's multiplicity, again according to Tien Chen, and reduces multiplicity 00885 by switching to the polynomial's derivative (which has the same root, with multiplicity 00886 reduced by one), as proposed by C. Bond.</li> 00887 </ul> 00888 00889 The algorithm has been successfully used for polynomials up to order 80. 00890 The function stops and returns <tt>false</tt> if an iteration fails to converge within 00891 80 steps. The type <tt>POLYNOMIAL</tt> must be compatible to 00892 \ref vigra::PolynomialView, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt> 00893 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>. 00894 00895 <b> Declaration:</b> 00896 00897 pass arguments explicitly: 00898 \code 00899 namespace vigra { 00900 template <class POLYNOMIAL, class VECTOR> 00901 bool 00902 polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots = true); 00903 } 00904 \endcode 00905 00906 00907 <b> Usage:</b> 00908 00909 <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br> 00910 Namespace: vigra 00911 00912 \code 00913 // encode the polynomial x^4 - 1 00914 Polynomial<double> poly(4); 00915 poly[0] = -1.0; 00916 poly[4] = 1.0; 00917 00918 ArrayVector<std::complex<double> > roots; 00919 polynomialRoots(poly, roots); 00920 \endcode 00921 00922 \see polynomialRootsEigenvalueMethod() 00923 */ 00924 template <class POLYNOMIAL, class VECTOR> 00925 bool polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots) 00926 { 00927 typedef typename POLYNOMIAL::value_type T; 00928 typedef typename POLYNOMIAL::Real Real; 00929 typedef typename POLYNOMIAL::Complex Complex; 00930 typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial; 00931 00932 double eps = poriginal.epsilon(); 00933 00934 WorkPolynomial p(poriginal.begin(), poriginal.order(), eps); 00935 p.minimizeOrder(); 00936 if(p.order() == 0) 00937 return true; 00938 00939 Complex x = detail::laguerreStartingGuess(p); 00940 00941 unsigned int multiplicity = 1; 00942 bool triedConjugate = false; 00943 00944 // handle the high order cases 00945 while(p.order() > 2) 00946 { 00947 p.balance(); 00948 00949 // find root estimate using Laguerre's method on deflated polynomial p; 00950 // zero return indicates failure to converge 00951 multiplicity = detail::laguerre1Root(p, x, multiplicity); 00952 00953 if(multiplicity == 0) 00954 return false; 00955 // polish root on original polynomial poriginal; 00956 // zero return indicates failure to converge 00957 if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity)) 00958 return false; 00959 x = detail::deleteBelowEpsilon(x, eps); 00960 roots.push_back(x); 00961 p.deflate(x); 00962 // determine the next starting guess 00963 if(multiplicity > 1) 00964 { 00965 // probably multiple root => keep current root as starting guess 00966 --multiplicity; 00967 triedConjugate = false; 00968 } 00969 else 00970 { 00971 // need a new starting guess 00972 if(x.imag() != 0.0 && !triedConjugate) 00973 { 00974 // if the root is complex and we don't already have 00975 // the conjugate root => try the conjugate as starting guess 00976 triedConjugate = true; 00977 x = conj(x); 00978 } 00979 else 00980 { 00981 // otherwise generate new starting guess 00982 triedConjugate = false; 00983 x = detail::laguerreStartingGuess(p); 00984 } 00985 } 00986 } 00987 00988 // handle the low order cases 00989 if(p.order() == 2) 00990 { 00991 Complex a = p[2]; 00992 Complex b = p[1]; 00993 Complex c = p[0]; 00994 Complex b2 = std::sqrt(b*b - 4.0*a*c); 00995 Complex q; 00996 if((conj(b)*b2).real() >= 0.0) 00997 q = -0.5 * (b + b2); 00998 else 00999 q = -0.5 * (b - b2); 01000 x = detail::complexDiv(q, a); 01001 if(polishRoots) 01002 detail::laguerre1Root(poriginal, x, 1); 01003 roots.push_back(detail::deleteBelowEpsilon(x, eps)); 01004 x = detail::complexDiv(c, q); 01005 if(polishRoots) 01006 detail::laguerre1Root(poriginal, x, 1); 01007 roots.push_back(detail::deleteBelowEpsilon(x, eps)); 01008 } 01009 else if(p.order() == 1) 01010 { 01011 x = detail::complexDiv(-p[0], p[1]); 01012 if(polishRoots) 01013 detail::laguerre1Root(poriginal, x, 1); 01014 roots.push_back(detail::deleteBelowEpsilon(x, eps)); 01015 } 01016 std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps)); 01017 return true; 01018 } 01019 01020 template <class POLYNOMIAL, class VECTOR> 01021 inline bool 01022 polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots) 01023 { 01024 return polynomialRoots(poriginal, roots, true); 01025 } 01026 01027 /** Determine the real roots of the polynomial <tt>p</tt>. 01028 01029 This function simply calls \ref polynomialRoots() and than throws away all complex roots. 01030 Accordingly, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt> 01031 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>. 01032 01033 <b> Declaration:</b> 01034 01035 pass arguments explicitly: 01036 \code 01037 namespace vigra { 01038 template <class POLYNOMIAL, class VECTOR> 01039 bool 01040 polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots = true); 01041 } 01042 \endcode 01043 01044 01045 <b> Usage:</b> 01046 01047 <b>\#include</b> "<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>"<br> 01048 Namespace: vigra 01049 01050 \code 01051 // encode the polynomial x^4 - 1 01052 Polynomial<double> poly(4); 01053 poly[0] = -1.0; 01054 poly[4] = 1.0; 01055 01056 ArrayVector<double> roots; 01057 polynomialRealRoots(poly, roots); 01058 \endcode 01059 01060 \see polynomialRealRootsEigenvalueMethod() 01061 */ 01062 template <class POLYNOMIAL, class VECTOR> 01063 bool polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots) 01064 { 01065 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex; 01066 ArrayVector<Complex> croots; 01067 if(!polynomialRoots(p, croots, polishRoots)) 01068 return false; 01069 for(unsigned int i = 0; i < croots.size(); ++i) 01070 if(croots[i].imag() == 0.0) 01071 roots.push_back(croots[i].real()); 01072 return true; 01073 } 01074 01075 template <class POLYNOMIAL, class VECTOR> 01076 inline bool 01077 polynomialRealRoots(POLYNOMIAL const & poriginal, VECTOR & roots) 01078 { 01079 return polynomialRealRoots(poriginal, roots, true); 01080 } 01081 01082 //@} 01083 01084 } // namespace vigra 01085 01086 #endif // VIGRA_POLYNOMIAL_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|