[ 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
49namespace vigra {
50
51template <class T> class Polynomial;
52template <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> <vigra/polynomial.hxx><br>
68 Namespace: vigra
69
70 \ingroup Polynomials
71*/
72template <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
102 typedef Polynomial<Real> RealPolynomial;
103 typedef Polynomial<Complex> ComplexPolynomial;
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 intermediate 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 smallest 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
263template <class T>
264template <class U>
265typename 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/*
277template <class T>
278typename PolynomialView<T>::Complex
279PolynomialView<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
290template <class T>
291void
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
310template <class T>
311void
312PolynomialView<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
327 }
328 if(multiplicity > 1)
329 deflate(v, multiplicity-1);
330}
331
332template <class T>
333void
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
344template <class T>
345void
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
365template <class T>
366void
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
378template <class T>
379void
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
396template <class T>
397void
399{
400 while(std::abs(coeffs_[order_]) <= epsilon && order_ > 0)
401 --order_;
402}
403
404template <class T>
405void
407{
408 for(unsigned int i = 0; i<order_; ++i)
409 coeffs_[i] /= coeffs_[order_];
410 coeffs_[order_] = T(1.0);
411}
412
413template <class T>
414void
415PolynomialView<T>::balance()
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> <vigra/polynomial.hxx><br>
438 Namespace: vigra
439
440 \ingroup Polynomials
441*/
442template <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;
450 typedef Polynomial<Real> RealPolynomial;
451 typedef Polynomial<Complex> ComplexPolynomial;
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> <vigra/polynomial.hxx><br>
569 Namespace: vigra
570
571 \ingroup Polynomials
572*/
573template <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;
582 typedef StaticPolynomial<MAXORDER, Real> RealPolynomial;
583 typedef StaticPolynomial<MAXORDER, Complex> ComplexPolynomial;
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
701namespace detail {
702
703// replacement for complex division (some compilers have numerically
704// less stable implementations); code from python complexobject.c
705template <class T>
706std::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
736template <class T>
737std::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
746template <class POLYNOMIAL>
747typename POLYNOMIAL::value_type
748laguerreStartingGuess(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
756template <class POLYNOMIAL, class Complex>
757int laguerre1Root(POLYNOMIAL const & p, Complex & x, unsigned int multiplicity)
758{
759 typedef typename NumericTraits<Complex>::ValueType Real;
760
761 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
858template <class Real>
859struct 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> <vigra/polynomial.hxx><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*/
938template <class POLYNOMIAL, class VECTOR>
939bool polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots, bool polishRoots)
940{
941 typedef typename POLYNOMIAL::Real Real;
942 typedef typename POLYNOMIAL::Complex Complex;
943 typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial;
944
945 double eps = poriginal.epsilon();
946
947 WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
948 p.minimizeOrder();
949 if(p.order() == 0)
950 return true;
951
952 Complex x = detail::laguerreStartingGuess(p);
953
954 unsigned int multiplicity = 1;
955 bool triedConjugate = false;
956
957 // handle the high order cases
958 while(p.order() > 2)
959 {
960 p.balance();
961
962 // find root estimate using Laguerre's method on deflated polynomial p;
963 // zero return indicates failure to converge
964 multiplicity = detail::laguerre1Root(p, x, multiplicity);
965
966 if(multiplicity == 0)
967 return false;
968 // polish root on original polynomial poriginal;
969 // zero return indicates failure to converge
970 if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
971 return false;
972 x = detail::deleteBelowEpsilon(x, eps);
973 roots.push_back(x);
974 p.deflate(x);
975 // determine the next starting guess
976 if(multiplicity > 1)
977 {
978 // probably multiple root => keep current root as starting guess
979 --multiplicity;
980 triedConjugate = false;
981 }
982 else
983 {
984 // need a new starting guess
985 if(x.imag() != 0.0 && !triedConjugate)
986 {
987 // if the root is complex and we don't already have
988 // the conjugate root => try the conjugate as starting guess
989 triedConjugate = true;
990 x = conj(x);
991 }
992 else
993 {
994 // otherwise generate new starting guess
995 triedConjugate = false;
996 x = detail::laguerreStartingGuess(p);
997 }
998 }
999 }
1000
1001 // handle the low order cases
1002 if(p.order() == 2)
1003 {
1004 Complex a = p[2];
1005 Complex b = p[1];
1006 Complex c = p[0];
1007 Complex b2 = std::sqrt(b*b - 4.0*a*c);
1008 Complex q;
1009 if((conj(b)*b2).real() >= 0.0)
1010 q = -0.5 * (b + b2);
1011 else
1012 q = -0.5 * (b - b2);
1013 x = detail::complexDiv(q, a);
1014 if(polishRoots)
1015 detail::laguerre1Root(poriginal, x, 1);
1016 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1017 x = detail::complexDiv(c, q);
1018 if(polishRoots)
1019 detail::laguerre1Root(poriginal, x, 1);
1020 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1021 }
1022 else if(p.order() == 1)
1023 {
1024 x = detail::complexDiv(-p[0], p[1]);
1025 if(polishRoots)
1026 detail::laguerre1Root(poriginal, x, 1);
1027 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1028 }
1029 std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
1030 return true;
1031}
1032
1033template <class POLYNOMIAL, class VECTOR>
1034inline bool
1035polynomialRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
1036{
1037 return polynomialRoots(poriginal, roots, true);
1038}
1039
1040/** Determine the real roots of the polynomial <tt>p</tt>.
1041
1042 This function simply calls \ref polynomialRoots() and than throws away all complex roots.
1043 Accordingly, <tt>VECTOR</tt> must be compatible to <tt>std::vector</tt>
1044 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>.
1045
1046 <b> Declaration:</b>
1047
1048 pass arguments explicitly:
1049 \code
1050 namespace vigra {
1051 template <class POLYNOMIAL, class VECTOR>
1052 bool
1053 polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots = true);
1054 }
1055 \endcode
1056
1057
1058 <b> Usage:</b>
1059
1060 <b>\#include</b> <vigra/polynomial.hxx><br>
1061 Namespace: vigra
1062
1063 \code
1064 // encode the polynomial x^4 - 1
1065 Polynomial<double> poly(4);
1066 poly[0] = -1.0;
1067 poly[4] = 1.0;
1068
1069 ArrayVector<double> roots;
1070 polynomialRealRoots(poly, roots);
1071 \endcode
1072
1073 \see polynomialRealRootsEigenvalueMethod()
1074*/
1075template <class POLYNOMIAL, class VECTOR>
1076bool polynomialRealRoots(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots)
1077{
1078 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1079 ArrayVector<Complex> croots;
1080 if(!polynomialRoots(p, croots, polishRoots))
1081 return false;
1082 for(unsigned int i = 0; i < croots.size(); ++i)
1083 if(croots[i].imag() == 0.0)
1084 roots.push_back(croots[i].real());
1085 return true;
1086}
1087
1088template <class POLYNOMIAL, class VECTOR>
1089inline bool
1090polynomialRealRoots(POLYNOMIAL const & poriginal, VECTOR & roots)
1091{
1092 return polynomialRealRoots(poriginal, roots, true);
1093}
1094
1095//@}
1096
1097} // namespace vigra
1098
1099namespace std {
1100
1101template <class T>
1102ostream & operator<<(ostream & o, vigra::PolynomialView<T> const & p)
1103{
1104 for(unsigned int k=0; k < p.order(); ++k)
1105 o << p[k] << " ";
1106 o << p[p.order()];
1107 return o;
1108}
1109
1110} // namespace std
1111
1112#endif // VIGRA_POLYNOMIAL_HXX
size_type size() const
Definition array_vector.hxx:358
Definition array_vector.hxx:514
void minimizeOrder(double epsilon=0.0)
Definition polynomial.hxx:398
PolynomialView(T *coeffs, unsigned int order, double epsilon=1.0e-14)
Definition polynomial.hxx:120
T const * const_iterator
Definition polynomial.hxx:100
T value_type
Definition polynomial.hxx:79
const_iterator begin() const
Definition polynomial.hxx:215
void backwardDeflate(T v)
Definition polynomial.hxx:367
NumericTraits< RealPromote >::ComplexPromote Complex
Definition polynomial.hxx:92
T * iterator
Definition polynomial.hxx:96
void deflate(T const &r, unsigned int multiplicity=1)
Definition polynomial.hxx:312
unsigned int size() const
Definition polynomial.hxx:225
double epsilon() const
Definition polynomial.hxx:236
void forwardDeflate(T const &v)
Definition polynomial.hxx:334
void deflateConjugatePair(Complex const &v)
Definition polynomial.hxx:380
T & operator[](unsigned int i)
Access the coefficient of x^i.
Definition polynomial.hxx:127
unsigned int order() const
Definition polynomial.hxx:230
void differentiate(unsigned int n=1)
Definition polynomial.hxx:292
void forwardBackwardDeflate(T v)
Definition polynomial.hxx:346
NumericTraits< RealPromote >::ValueType Real
Definition polynomial.hxx:88
void setEpsilon(double eps)
Definition polynomial.hxx:242
NumericTraits< T >::RealPromote RealPromote
Definition polynomial.hxx:84
iterator end()
Definition polynomial.hxx:210
const_iterator end() const
Definition polynomial.hxx:220
void normalize()
Definition polynomial.hxx:406
iterator begin()
Definition polynomial.hxx:205
T const & operator[](unsigned int i) const
Access the coefficient of x^i.
Definition polynomial.hxx:131
PromoteTraits< T, V >::Promote operator()(V v) const
Definition polynomial.hxx:445
Polynomial(Polynomial const &p)
Definition polynomial.hxx:472
Polynomial< Complex > getDeflated(Complex const &r) const
Definition polynomial.hxx:542
Polynomial(unsigned int order=0, double epsilon=1.0e-14)
Definition polynomial.hxx:463
Polynomial< T > getDerivative(unsigned int n=1) const
Definition polynomial.hxx:518
Polynomial & operator=(Polynomial const &p)
Definition polynomial.hxx:504
Polynomial< T > getDeflated(Real r) const
Definition polynomial.hxx:529
Polynomial(ITER i, unsigned int order, double epsilon)
Definition polynomial.hxx:495
Polynomial(ITER i, unsigned int order)
Definition polynomial.hxx:482
Definition polynomial.hxx:576
StaticPolynomial getDeflated(Real r) const
Definition polynomial.hxx:668
StaticPolynomial(StaticPolynomial const &p)
Definition polynomial.hxx:607
StaticPolynomial getDerivative(unsigned int n=1) const
Definition polynomial.hxx:657
StaticPolynomial< MAXORDER, Complex > getDeflated(Complex const &r) const
Definition polynomial.hxx:681
StaticPolynomial & operator=(StaticPolynomial const &p)
Definition polynomial.hxx:644
StaticPolynomial(unsigned int order=0, double epsilon=1.0e-14)
Definition polynomial.hxx:596
StaticPolynomial(ITER i, unsigned int order)
Definition polynomial.hxx:618
StaticPolynomial(ITER i, unsigned int order, double epsilon)
Definition polynomial.hxx:633
R imag(const FFTWComplex< R > &a)
imaginary part
Definition fftw3.hxx:1023
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition fftw3.hxx:1030
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
R real(const FFTWComplex< R > &a)
real part
Definition fftw3.hxx:1016
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition fftw3.hxx:1002
bool polynomialRoots(POLYNOMIAL const &poriginal, VECTOR &roots, bool polishRoots)
Definition polynomial.hxx:939
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point equality.
Definition mathutil.hxx:1638
FixedPoint< 0, FracBits > frac(FixedPoint< IntBits, FracBits > v)
fractional part.
Definition fixedpoint.hxx:650
bool polynomialRealRoots(POLYNOMIAL const &p, VECTOR &roots, bool polishRoots)
Definition polynomial.hxx:1076

© 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.12.1 (Thu Feb 27 2025)