[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

polynomial.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_POLYNOMIAL_HXX
38 #define VIGRA_POLYNOMIAL_HXX
39 
40 #include <cmath>
41 #include <complex>
42 #include <algorithm>
43 #include <iosfwd>
44 #include "error.hxx"
45 #include "mathutil.hxx"
46 #include "numerictraits.hxx"
47 #include "array_vector.hxx"
48 
49 namespace vigra {
50 
51 template <class T> class Polynomial;
52 template <unsigned int MAXORDER, class T> class StaticPolynomial;
53 
54 /*****************************************************************/
55 /* */
56 /* PolynomialView */
57 /* */
58 /*****************************************************************/
59 
60 /** Polynomial interface for an externally managed array.
61 
62  The coefficient type <tt>T</tt> can be either a scalar or complex
63  (compatible to <tt>std::complex</tt>) type.
64 
65  \see vigra::Polynomial, vigra::StaticPolynomial, polynomialRoots()
66 
67  <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
68  Namespace: vigra
69 
70  \ingroup Polynomials
71 */
72 template <class T>
74 {
75  public:
76 
77  /** Coefficient type of the polynomial
78  */
79  typedef T value_type;
80 
81  /** Promote type of <tt>value_type</tt>
82  used for polynomial calculations
83  */
84  typedef typename NumericTraits<T>::RealPromote RealPromote;
85 
86  /** Scalar type associated with <tt>RealPromote</tt>
87  */
88  typedef typename NumericTraits<RealPromote>::ValueType Real;
89 
90  /** Complex type associated with <tt>RealPromote</tt>
91  */
92  typedef typename NumericTraits<RealPromote>::ComplexPromote Complex;
93 
94  /** Iterator for the coefficient sequence
95  */
96  typedef T * iterator;
97 
98  /** Const iterator for the coefficient sequence
99  */
100  typedef T const * const_iterator;
101 
104 
105 
106  /** Construct from a coefficient array of given <tt>order</tt>.
107 
108  The externally managed array must have length <tt>order+1</tt>
109  and is interpreted as representing the polynomial:
110 
111  \code
112  coeffs[0] + x * coeffs[1] + x * x * coeffs[2] + ...
113  \endcode
114 
115  The coefficients are not copied, we only store a pointer to the
116  array.<tt>epsilon</tt> (default: 1.0e-14) determines the precision
117  of subsequent algorithms (especially root finding) performed on the
118  polynomial.
119  */
120  PolynomialView(T * coeffs, unsigned int order, double epsilon = 1.0e-14)
121  : coeffs_(coeffs),
122  order_(order),
123  epsilon_(epsilon)
124  {}
125 
126  /// Access the coefficient of x^i
127  T & operator[](unsigned int i)
128  { return coeffs_[i]; }
129 
130  /// Access the coefficient of x^i
131  T const & operator[](unsigned int i) const
132  { return coeffs_[i]; }
133 
134  /** Evaluate the polynomial at the point <tt>v</tt>
135 
136  Multiplication must be defined between the types
137  <tt>V</tt> and <tt>PromoteTraits<T, V>::Promote</tt>.
138  If both <tt>V</tt> and <tt>V</tt> are scalar, the result will
139  be a scalar, otherwise it will be complex.
140  */
141  template <class V>
142  typename PromoteTraits<T, V>::Promote
143  operator()(V v) const;
144 
145  /** Differentiate the polynomial <tt>n</tt> times.
146  */
147  void differentiate(unsigned int n = 1);
148 
149  /** Deflate the polynomial at the root <tt>r</tt> with
150  the given <tt>multiplicity</tt>.
151 
152  The behavior of this function is undefined if <tt>r</tt>
153  is not a root with at least the given multiplicity.
154  This function calls forwardBackwardDeflate().
155  */
156  void deflate(T const & r, unsigned int multiplicity = 1);
157 
158  /** Forward-deflate the polynomial at the root <tt>r</tt>.
159 
160  The behavior of this function is undefined if <tt>r</tt>
161  is not a root. Forward deflation is best if <tt>r</tt> is
162  the biggest root (by magnitude).
163  */
164  void forwardDeflate(T const & v);
165 
166  /** Forward/backward eflate the polynomial at the root <tt>r</tt>.
167 
168  The behavior of this function is undefined if <tt>r</tt>
169  is not a root. Combined forward/backward deflation is best
170  if <tt>r</tt> is an ontermediate root or we don't know
171  <tt>r</tt>'s relation to the other roots of the polynomial.
172  */
173  void forwardBackwardDeflate(T v);
174 
175  /** Backward-deflate the polynomial at the root <tt>r</tt>.
176 
177  The behavior of this function is undefined if <tt>r</tt>
178  is not a root. Backward deflation is best if <tt>r</tt> is
179  the snallest root (by magnitude).
180  */
181  void backwardDeflate(T v);
182 
183  /** Deflate the polynomial with the complex conjugate roots
184  <tt>r</tt> and <tt>conj(r)</tt>.
185 
186  The behavior of this function is undefined if these are not
187  roots.
188  */
189  void deflateConjugatePair(Complex const & v);
190 
191  /** Adjust the polynomial's order if the highest coefficients are near zero.
192  The order is reduced as long as the absolute value does not exceed
193  the given \a epsilon.
194  */
195  void minimizeOrder(double epsilon = 0.0);
196 
197  /** Normalize the polynomial, i.e. dived by the highest coefficient.
198  */
199  void normalize();
200 
201  void balance();
202 
203  /** Get iterator for the coefficient sequence.
204  */
206  { return coeffs_; }
207 
208  /** Get end iterator for the coefficient sequence.
209  */
211  { return begin() + size(); }
212 
213  /** Get const_iterator for the coefficient sequence.
214  */
216  { return coeffs_; }
217 
218  /** Get end const_iterator for the coefficient sequence.
219  */
221  { return begin() + size(); }
222 
223  /** Get length of the coefficient sequence (<tt>order() + 1</tt>).
224  */
225  unsigned int size() const
226  { return order_ + 1; }
227 
228  /** Get order of the polynomial.
229  */
230  unsigned int order() const
231  { return order_; }
232 
233  /** Get requested precision for polynomial algorithms
234  (especially root finding).
235  */
236  double epsilon() const
237  { return epsilon_; }
238 
239  /** Set requested precision for polynomial algorithms
240  (especially root finding).
241  */
242  void setEpsilon(double eps)
243  { epsilon_ = eps; }
244 
245  protected:
246  PolynomialView(double epsilon = 1e-14)
247  : coeffs_(0),
248  order_(0),
249  epsilon_(epsilon)
250  {}
251 
252  void setCoeffs(T * coeffs, unsigned int order)
253  {
254  coeffs_ = coeffs;
255  order_ = order;
256  }
257 
258  T * coeffs_;
259  unsigned int order_;
260  double epsilon_;
261 };
262 
263 template <class T>
264 template <class U>
265 typename PromoteTraits<T, U>::Promote
267 {
268  typename PromoteTraits<T, U>::Promote p(coeffs_[order_]);
269  for(int i = order_ - 1; i >= 0; --i)
270  {
271  p = v * p + coeffs_[i];
272  }
273  return p;
274 }
275 
276 /*
277 template <class T>
278 typename PolynomialView<T>::Complex
279 PolynomialView<T>::operator()(Complex const & v) const
280 {
281  Complex p(coeffs_[order_]);
282  for(int i = order_ - 1; i >= 0; --i)
283  {
284  p = v * p + coeffs_[i];
285  }
286  return p;
287 }
288 */
289 
290 template <class T>
291 void
293 {
294  if(n == 0)
295  return;
296  if(order_ == 0)
297  {
298  coeffs_[0] = 0.0;
299  return;
300  }
301  for(unsigned int i = 1; i <= order_; ++i)
302  {
303  coeffs_[i-1] = double(i)*coeffs_[i];
304  }
305  --order_;
306  if(n > 1)
307  differentiate(n-1);
308 }
309 
310 template <class T>
311 void
312 PolynomialView<T>::deflate(T const & v, unsigned int multiplicity)
313 {
314  vigra_precondition(order_ > 0,
315  "PolynomialView<T>::deflate(): cannot deflate 0th order polynomial.");
316  if(v == 0.0)
317  {
318  ++coeffs_;
319  --order_;
320  }
321  else
322  {
323  // we use combined forward/backward deflation because
324  // our initial guess seems to favour convergence to
325  // a root with magnitude near the median among all roots
326  forwardBackwardDeflate(v);
327  }
328  if(multiplicity > 1)
329  deflate(v, multiplicity-1);
330 }
331 
332 template <class T>
333 void
335 {
336  for(int i = order_-1; i > 0; --i)
337  {
338  coeffs_[i] += v * coeffs_[i+1];
339  }
340  ++coeffs_;
341  --order_;
342 }
343 
344 template <class T>
345 void
347 {
348  unsigned int order2 = order_ / 2;
349  T tmp = coeffs_[order_];
350  for(unsigned int i = order_-1; i >= order2; --i)
351  {
352  T tmp1 = coeffs_[i];
353  coeffs_[i] = tmp;
354  tmp = tmp1 + v * tmp;
355  }
356  v = -1.0 / v;
357  coeffs_[0] *= v;
358  for(unsigned int i = 1; i < order2; ++i)
359  {
360  coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
361  }
362  --order_;
363 }
364 
365 template <class T>
366 void
368 {
369  v = -1.0 / v;
370  coeffs_[0] *= v;
371  for(unsigned int i = 1; i < order_; ++i)
372  {
373  coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
374  }
375  --order_;
376 }
377 
378 template <class T>
379 void
381 {
382  vigra_precondition(order_ > 1,
383  "PolynomialView<T>::deflateConjugatePair(): cannot deflate 2 roots "
384  "from 1st order polynomial.");
385  Real a = 2.0*v.real();
386  Real b = -sq(v.real()) - sq(v.imag());
387  coeffs_[order_-1] += a * coeffs_[order_];
388  for(int i = order_-2; i > 1; --i)
389  {
390  coeffs_[i] += a * coeffs_[i+1] + b*coeffs_[i+2];
391  }
392  coeffs_ += 2;
393  order_ -= 2;
394 }
395 
396 template <class T>
397 void
399 {
400  while(std::abs(coeffs_[order_]) <= epsilon && order_ > 0)
401  --order_;
402 }
403 
404 template <class T>
405 void
407 {
408  for(unsigned int i = 0; i<order_; ++i)
409  coeffs_[i] /= coeffs_[order_];
410  coeffs_[order_] = T(1.0);
411 }
412 
413 template <class T>
414 void
416 {
417  Real p0 = abs(coeffs_[0]), po = abs(coeffs_[order_]);
418  Real norm = (p0 > 0.0)
419  ? VIGRA_CSTD::sqrt(p0*po)
420  : po;
421  for(unsigned int i = 0; i<=order_; ++i)
422  coeffs_[i] /= norm;
423 }
424 
425 /*****************************************************************/
426 /* */
427 /* Polynomial */
428 /* */
429 /*****************************************************************/
430 
431 /** Polynomial with internally managed array.
432 
433  Most interesting functionality is inherited from \ref vigra::PolynomialView.
434 
435  \see vigra::PolynomialView, vigra::StaticPolynomial, polynomialRoots()
436 
437  <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
438  Namespace: vigra
439 
440  \ingroup Polynomials
441 */
442 template <class T>
444 : public PolynomialView<T>
445 {
446  typedef PolynomialView<T> BaseType;
447  public:
448  typedef typename BaseType::Real Real;
449  typedef typename BaseType::Complex Complex;
452 
453  typedef T value_type;
454  typedef T * iterator;
455  typedef T const * const_iterator;
456 
457  /** Construct polynomial with given <tt>order</tt> and all coefficients
458  set to zero (they can be set later using <tt>operator[]</tt>
459  or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines
460  the precision of subsequent algorithms (especially root finding)
461  performed on the polynomial.
462  */
463  Polynomial(unsigned int order = 0, double epsilon = 1.0e-14)
464  : BaseType(epsilon),
465  polynomial_(order + 1, T())
466  {
467  this->setCoeffs(&polynomial_[0], order);
468  }
469 
470  /** Copy constructor
471  */
473  : BaseType(p.epsilon()),
474  polynomial_(p.begin(), p.end())
475  {
476  this->setCoeffs(&polynomial_[0], p.order());
477  }
478 
479  /** Construct polynomial by copying the given coefficient sequence.
480  */
481  template <class ITER>
482  Polynomial(ITER i, unsigned int order)
483  : BaseType(),
484  polynomial_(i, i + order + 1)
485  {
486  this->setCoeffs(&polynomial_[0], order);
487  }
488 
489  /** Construct polynomial by copying the given coefficient sequence.
490  Set <tt>epsilon</tt> (default: 1.0e-14) as
491  the precision of subsequent algorithms (especially root finding)
492  performed on the polynomial.
493  */
494  template <class ITER>
495  Polynomial(ITER i, unsigned int order, double epsilon)
496  : BaseType(epsilon),
497  polynomial_(i, i + order + 1)
498  {
499  this->setCoeffs(&polynomial_[0], order);
500  }
501 
502  /** Assigment
503  */
505  {
506  if(this == &p)
507  return *this;
508  ArrayVector<T> tmp(p.begin(), p.end());
509  polynomial_.swap(tmp);
510  this->setCoeffs(&polynomial_[0], p.order());
511  this->epsilon_ = p.epsilon_;
512  return *this;
513  }
514 
515  /** Construct new polynomial representing the derivative of this
516  polynomial.
517  */
518  Polynomial<T> getDerivative(unsigned int n = 1) const
519  {
520  Polynomial<T> res(*this);
521  res.differentiate(n);
522  return res;
523  }
524 
525  /** Construct new polynomial representing this polynomial after
526  deflation at the real root <tt>r</tt>.
527  */
529  getDeflated(Real r) const
530  {
531  Polynomial<T> res(*this);
532  res.deflate(r);
533  return res;
534  }
535 
536  /** Construct new polynomial representing this polynomial after
537  deflation at the complex root <tt>r</tt>. The resulting
538  polynomial will have complex coefficients, even if this
539  polynomial had real ones.
540  */
542  getDeflated(Complex const & r) const
543  {
544  Polynomial<Complex> res(this->begin(), this->order(), this->epsilon());
545  res.deflate(r);
546  return res;
547  }
548 
549  protected:
550  ArrayVector<T> polynomial_;
551 };
552 
553 /*****************************************************************/
554 /* */
555 /* StaticPolynomial */
556 /* */
557 /*****************************************************************/
558 
559 /** Polynomial with internally managed array of static length.
560 
561  Most interesting functionality is inherited from \ref vigra::PolynomialView.
562  This class differs from \ref vigra::Polynomial in that it allocates
563  its memory statically which is much faster. Therefore, <tt>StaticPolynomial</tt>
564  can only represent polynomials up to the given <tt>MAXORDER</tt>.
565 
566  \see vigra::PolynomialView, vigra::Polynomial, polynomialRoots()
567 
568  <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
569  Namespace: vigra
570 
571  \ingroup Polynomials
572 */
573 template <unsigned int MAXORDER, class T>
575 : public PolynomialView<T>
576 {
577  typedef PolynomialView<T> BaseType;
578 
579  public:
580  typedef typename BaseType::Real Real;
581  typedef typename BaseType::Complex Complex;
584 
585  typedef T value_type;
586  typedef T * iterator;
587  typedef T const * const_iterator;
588 
589 
590  /** Construct polynomial with given <tt>order <= MAXORDER</tt> and all
591  coefficients set to zero (they can be set later using <tt>operator[]</tt>
592  or the iterators). <tt>epsilon</tt> (default: 1.0e-14) determines
593  the precision of subsequent algorithms (especially root finding)
594  performed on the polynomial.
595  */
596  StaticPolynomial(unsigned int order = 0, double epsilon = 1.0e-14)
597  : BaseType(epsilon)
598  {
599  vigra_precondition(order <= MAXORDER,
600  "StaticPolynomial(): order exceeds MAXORDER.");
601  std::fill_n(polynomial_, order+1, T());
602  this->setCoeffs(polynomial_, order);
603  }
604 
605  /** Copy constructor
606  */
608  : BaseType(p.epsilon())
609  {
610  std::copy(p.begin(), p.end(), polynomial_);
611  this->setCoeffs(polynomial_, p.order());
612  }
613 
614  /** Construct polynomial by copying the given coefficient sequence.
615  <tt>order <= MAXORDER</tt> is required.
616  */
617  template <class ITER>
618  StaticPolynomial(ITER i, unsigned int order)
619  : BaseType()
620  {
621  vigra_precondition(order <= MAXORDER,
622  "StaticPolynomial(): order exceeds MAXORDER.");
623  std::copy(i, i + order + 1, polynomial_);
624  this->setCoeffs(polynomial_, order);
625  }
626 
627  /** Construct polynomial by copying the given coefficient sequence.
628  <tt>order <= MAXORDER</tt> is required. Set <tt>epsilon</tt> (default: 1.0e-14) as
629  the precision of subsequent algorithms (especially root finding)
630  performed on the polynomial.
631  */
632  template <class ITER>
633  StaticPolynomial(ITER i, unsigned int order, double epsilon)
634  : BaseType(epsilon)
635  {
636  vigra_precondition(order <= MAXORDER,
637  "StaticPolynomial(): order exceeds MAXORDER.");
638  std::copy(i, i + order + 1, polynomial_);
639  this->setCoeffs(polynomial_, order);
640  }
641 
642  /** Assigment.
643  */
645  {
646  if(this == &p)
647  return *this;
648  std::copy(p.begin(), p.end(), polynomial_);
649  this->setCoeffs(polynomial_, p.order());
650  this->epsilon_ = p.epsilon_;
651  return *this;
652  }
653 
654  /** Construct new polynomial representing the derivative of this
655  polynomial.
656  */
657  StaticPolynomial getDerivative(unsigned int n = 1) const
658  {
659  StaticPolynomial res(*this);
660  res.differentiate(n);
661  return res;
662  }
663 
664  /** Construct new polynomial representing this polynomial after
665  deflation at the real root <tt>r</tt>.
666  */
668  getDeflated(Real r) const
669  {
670  StaticPolynomial res(*this);
671  res.deflate(r);
672  return res;
673  }
674 
675  /** Construct new polynomial representing this polynomial after
676  deflation at the complex root <tt>r</tt>. The resulting
677  polynomial will have complex coefficients, even if this
678  polynomial had real ones.
679  */
681  getDeflated(Complex const & r) const
682  {
683  StaticPolynomial<MAXORDER, Complex> res(this->begin(), this->order(), this->epsilon());
684  res.deflate(r);
685  return res;
686  }
687 
688  void setOrder(unsigned int order)
689  {
690  vigra_precondition(order <= MAXORDER,
691  "taticPolynomial::setOrder(): order exceeds MAXORDER.");
692  this->order_ = order;
693  }
694 
695  protected:
696  T polynomial_[MAXORDER+1];
697 };
698 
699 /************************************************************/
700 
701 namespace detail {
702 
703 // replacement for complex division (some compilers have numerically
704 // less stable implementations); code from python complexobject.c
705 template <class T>
706 std::complex<T> complexDiv(std::complex<T> const & a, std::complex<T> const & b)
707 {
708  const double abs_breal = b.real() < 0 ? -b.real() : b.real();
709  const double abs_bimag = b.imag() < 0 ? -b.imag() : b.imag();
710 
711  if (abs_breal >= abs_bimag)
712  {
713  /* divide tops and bottom by b.real() */
714  if (abs_breal == 0.0)
715  {
716  return std::complex<T>(a.real() / abs_breal, a.imag() / abs_breal);
717  }
718  else
719  {
720  const double ratio = b.imag() / b.real();
721  const double denom = b.real() + b.imag() * ratio;
722  return std::complex<T>((a.real() + a.imag() * ratio) / denom,
723  (a.imag() - a.real() * ratio) / denom);
724  }
725  }
726  else
727  {
728  /* divide tops and bottom by b.imag() */
729  const double ratio = b.real() / b.imag();
730  const double denom = b.real() * ratio + b.imag();
731  return std::complex<T>((a.real() * ratio + a.imag()) / denom,
732  (a.imag() * ratio - a.real()) / denom);
733  }
734 }
735 
736 template <class T>
737 std::complex<T> deleteBelowEpsilon(std::complex<T> const & x, double eps)
738 {
739  return std::abs(x.imag()) <= 2.0*eps*std::abs(x.real())
740  ? std::complex<T>(x.real())
741  : std::abs(x.real()) <= 2.0*eps*std::abs(x.imag())
742  ? std::complex<T>(NumericTraits<T>::zero(), x.imag())
743  : x;
744 }
745 
746 template <class POLYNOMIAL>
747 typename POLYNOMIAL::value_type
748 laguerreStartingGuess(POLYNOMIAL const & p)
749 {
750  double N = p.order();
751  typename POLYNOMIAL::value_type centroid = -p[p.order()-1] / N / p[p.order()];
752  double dist = VIGRA_CSTD::pow(std::abs(p(centroid) / p[p.order()]), 1.0 / N);
753  return centroid + dist;
754 }
755 
756 template <class POLYNOMIAL, class Complex>
757 int laguerre1Root(POLYNOMIAL const & p, Complex & x, unsigned int multiplicity)
758 {
759  typedef typename NumericTraits<Complex>::ValueType Real;
760 
761  static double frac[] = {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
762  int maxiter = 80,
763  count;
764  double N = p.order();
765  double eps = p.epsilon(),
766  eps2 = VIGRA_CSTD::sqrt(eps);
767 
768  if(multiplicity == 0)
769  x = laguerreStartingGuess(p);
770 
771  bool mayTryDerivative = true; // try derivative for multiple roots
772 
773  for(count = 0; count < maxiter; ++count)
774  {
775  // Horner's algorithm to calculate values of polynomial and its
776  // first two derivatives and estimate error for current x
777  Complex p0(p[p.order()]);
778  Complex p1(0.0);
779  Complex p2(0.0);
780  Real ax = std::abs(x);
781  Real err = std::abs(p0);
782  for(int i = p.order()-1; i >= 0; --i)
783  {
784  p2 = p2 * x + p1;
785  p1 = p1 * x + p0;
786  p0 = p0 * x + p[i];
787  err = err * ax + std::abs(p0);
788  }
789  p2 *= 2.0;
790  err *= eps;
791  Real ap0 = std::abs(p0);
792  if(ap0 <= err)
793  {
794  break; // converged
795  }
796  Complex g = complexDiv(p1, p0);
797  Complex g2 = g * g;
798  Complex h = g2 - complexDiv(p2, p0);
799  // estimate root multiplicity according to Tien Chen
800  if(g2 != 0.0)
801  {
802  multiplicity = (unsigned int)VIGRA_CSTD::floor(N /
803  (std::abs(N * complexDiv(h, g2) - 1.0) + 1.0) + 0.5);
804  if(multiplicity < 1)
805  multiplicity = 1;
806  }
807  // improve accuracy of multiple roots on the derivative, as suggested by C. Bond
808  // (do this only if we are already near the root, otherwise we may converge to
809  // a different root of the derivative polynomial)
810  if(mayTryDerivative && multiplicity > 1 && ap0 < eps2)
811  {
812  Complex x1 = x;
813  int derivativeMultiplicity = laguerre1Root(p.getDerivative(), x1, multiplicity-1);
814  if(derivativeMultiplicity && std::abs(p(x1)) < std::abs(p(x)))
815  {
816  // successful search on derivative
817  x = x1;
818  return derivativeMultiplicity + 1;
819  }
820  else
821  {
822  // unsuccessful search on derivative => don't do it again
823  mayTryDerivative = false;
824  }
825  }
826  Complex sq = VIGRA_CSTD::sqrt((N - 1.0) * (N * h - g2));
827  Complex gp = g + sq;
828  Complex gm = g - sq;
829  if(std::abs(gp) < std::abs(gm))
830  gp = gm;
831  Complex dx;
832  if(gp != 0.0)
833  {
834  dx = complexDiv(Complex(N) , gp);
835  }
836  else
837  {
838  // re-initialisation trick due to Numerical Recipes
839  dx = (1.0 + ax) * Complex(VIGRA_CSTD::cos(double(count)), VIGRA_CSTD::sin(double(count)));
840  }
841  Complex x1 = x - dx;
842 
843  if(x1 - x == 0.0)
844  {
845  break; // converged
846  }
847  if((count + 1) % 10)
848  x = x1;
849  else
850  // cycle breaking trick according to Numerical Recipes
851  x = x - frac[(count+1)/10] * dx;
852  }
853  return count < maxiter ?
854  multiplicity :
855  0;
856 }
857 
858 template <class Real>
859 struct PolynomialRootCompare
860 {
861  Real epsilon;
862 
863  PolynomialRootCompare(Real eps)
864  : epsilon(eps)
865  {}
866 
867  template <class T>
868  bool operator()(T const & l, T const & r)
869  {
870  return closeAtTolerance(l.real(), r.real(), epsilon)
871  ? l.imag() < r.imag()
872  : l.real() < r.real();
873  }
874 };
875 
876 } // namespace detail
877 
878 /** \addtogroup Polynomials Polynomials and root determination
879 
880  Classes to represent polynomials and functions to find polynomial roots.
881 */
882 //@{
883 
884 /*****************************************************************/
885 /* */
886 /* polynomialRoots */
887 /* */
888 /*****************************************************************/
889 
890 /** Determine the roots of the polynomial <tt>poriginal</tt>.
891 
892  The roots are appended to the vector <tt>roots</tt>, with optional root
893  polishing as specified by <tt>polishRoots</tt> (default: do polishing). The function uses an
894  improved version of Laguerre's algorithm. The improvements are as follows:
895 
896  <ul>
897  <li>It uses a clever initial guess for the iteration, according to a proposal by Tien Chen</li>
898  <li>It estimates each root's multiplicity, again according to Tien Chen, and reduces multiplicity
899  by switching to the polynomial's derivative (which has the same root, with multiplicity
900  reduced by one), as proposed by C. Bond.</li>
901  </ul>
902 
903  The algorithm has been successfully used for polynomials up to order 80.
904  The function stops and returns <tt>false</tt> if an iteration fails to converge within
905  80 steps. The type <tt>POLYNOMIAL</tt> must be compatible to
906  \ref vigra::PolynomialView, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
907  with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>.
908 
909  <b> Declaration:</b>
910 
911  pass arguments explicitly:
912  \code
913  namespace vigra {
914  template <class POLYNOMIAL, class VECTOR>
915  bool
916  polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots = true);
917  }
918  \endcode
919 
920 
921  <b> Usage:</b>
922 
923  <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
924  Namespace: vigra
925 
926  \code
927  // encode the polynomial x^4 - 1
928  Polynomial<double> poly(4);
929  poly[0] = -1.0;
930  poly[4] = 1.0;
931 
932  ArrayVector<std::complex<double> > roots;
933  polynomialRoots(poly, roots);
934  \endcode
935 
936  \see polynomialRootsEigenvalueMethod()
937 */
938 template <class POLYNOMIAL, class VECTOR>
939 bool polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots)
940 {
941  typedef typename POLYNOMIAL::value_type T;
942  typedef typename POLYNOMIAL::Real Real;
943  typedef typename POLYNOMIAL::Complex Complex;
944  typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial;
945 
946  double eps = poriginal.epsilon();
947 
948  WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
949  p.minimizeOrder();
950  if(p.order() == 0)
951  return true;
952 
953  Complex x = detail::laguerreStartingGuess(p);
954 
955  unsigned int multiplicity = 1;
956  bool triedConjugate = false;
957 
958  // handle the high order cases
959  while(p.order() > 2)
960  {
961  p.balance();
962 
963  // find root estimate using Laguerre's method on deflated polynomial p;
964  // zero return indicates failure to converge
965  multiplicity = detail::laguerre1Root(p, x, multiplicity);
966 
967  if(multiplicity == 0)
968  return false;
969  // polish root on original polynomial poriginal;
970  // zero return indicates failure to converge
971  if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
972  return false;
973  x = detail::deleteBelowEpsilon(x, eps);
974  roots.push_back(x);
975  p.deflate(x);
976  // determine the next starting guess
977  if(multiplicity > 1)
978  {
979  // probably multiple root => keep current root as starting guess
980  --multiplicity;
981  triedConjugate = false;
982  }
983  else
984  {
985  // need a new starting guess
986  if(x.imag() != 0.0 && !triedConjugate)
987  {
988  // if the root is complex and we don't already have
989  // the conjugate root => try the conjugate as starting guess
990  triedConjugate = true;
991  x = conj(x);
992  }
993  else
994  {
995  // otherwise generate new starting guess
996  triedConjugate = false;
997  x = detail::laguerreStartingGuess(p);
998  }
999  }
1000  }
1001 
1002  // handle the low order cases
1003  if(p.order() == 2)
1004  {
1005  Complex a = p[2];
1006  Complex b = p[1];
1007  Complex c = p[0];
1008  Complex b2 = std::sqrt(b*b - 4.0*a*c);
1009  Complex q;
1010  if((conj(b)*b2).real() >= 0.0)
1011  q = -0.5 * (b + b2);
1012  else
1013  q = -0.5 * (b - b2);
1014  x = detail::complexDiv(q, a);
1015  if(polishRoots)
1016  detail::laguerre1Root(poriginal, x, 1);
1017  roots.push_back(detail::deleteBelowEpsilon(x, eps));
1018  x = detail::complexDiv(c, q);
1019  if(polishRoots)
1020  detail::laguerre1Root(poriginal, x, 1);
1021  roots.push_back(detail::deleteBelowEpsilon(x, eps));
1022  }
1023  else if(p.order() == 1)
1024  {
1025  x = detail::complexDiv(-p[0], p[1]);
1026  if(polishRoots)
1027  detail::laguerre1Root(poriginal, x, 1);
1028  roots.push_back(detail::deleteBelowEpsilon(x, eps));
1029  }
1030  std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
1031  return true;
1032 }
1033 
1034 template <class POLYNOMIAL, class VECTOR>
1035 inline bool
1036 polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
1037 {
1038  return polynomialRoots(poriginal, roots, true);
1039 }
1040 
1041 /** Determine the real roots of the polynomial <tt>p</tt>.
1042 
1043  This function simply calls \ref polynomialRoots() and than throws away all complex roots.
1044  Accordingly, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
1045  with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>.
1046 
1047  <b> Declaration:</b>
1048 
1049  pass arguments explicitly:
1050  \code
1051  namespace vigra {
1052  template <class POLYNOMIAL, class VECTOR>
1053  bool
1054  polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots = true);
1055  }
1056  \endcode
1057 
1058 
1059  <b> Usage:</b>
1060 
1061  <b>\#include</b> <<a href="polynomial_8hxx-source.html">vigra/polynomial.hxx</a>><br>
1062  Namespace: vigra
1063 
1064  \code
1065  // encode the polynomial x^4 - 1
1066  Polynomial<double> poly(4);
1067  poly[0] = -1.0;
1068  poly[4] = 1.0;
1069 
1070  ArrayVector<double> roots;
1071  polynomialRealRoots(poly, roots);
1072  \endcode
1073 
1074  \see polynomialRealRootsEigenvalueMethod()
1075 */
1076 template <class POLYNOMIAL, class VECTOR>
1077 bool polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots)
1078 {
1079  typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1080  ArrayVector<Complex> croots;
1081  if(!polynomialRoots(p, croots, polishRoots))
1082  return false;
1083  for(unsigned int i = 0; i < croots.size(); ++i)
1084  if(croots[i].imag() == 0.0)
1085  roots.push_back(croots[i].real());
1086  return true;
1087 }
1088 
1089 template <class POLYNOMIAL, class VECTOR>
1090 inline bool
1091 polynomialRealRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
1092 {
1093  return polynomialRealRoots(poriginal, roots, true);
1094 }
1095 
1096 //@}
1097 
1098 } // namespace vigra
1099 
1100 namespace std {
1101 
1102 template <class T>
1103 ostream & operator<<(ostream & o, vigra::PolynomialView<T> const & p)
1104 {
1105  for(unsigned int k=0; k < p.order(); ++k)
1106  o << p[k] << " ";
1107  o << p[p.order()];
1108  return o;
1109 }
1110 
1111 } // namespace std
1112 
1113 #endif // VIGRA_POLYNOMIAL_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.1 (Thu Jun 14 2012)