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