37 #ifndef VIGRA_POLYNOMIAL_HXX
38 #define VIGRA_POLYNOMIAL_HXX
45 #include "mathutil.hxx"
46 #include "numerictraits.hxx"
47 #include "array_vector.hxx"
51 template <
class T>
class Polynomial;
52 template <
unsigned int MAXORDER,
class T>
class StaticPolynomial;
88 typedef typename NumericTraits<RealPromote>::ValueType
Real;
92 typedef typename NumericTraits<RealPromote>::ComplexPromote
Complex;
128 {
return coeffs_[i]; }
132 {
return coeffs_[i]; }
142 typename PromoteTraits<T, V>::Promote
156 void deflate(T
const & r,
unsigned int multiplicity = 1);
226 {
return order_ + 1; }
252 void setCoeffs(T * coeffs,
unsigned int order)
265 typename PromoteTraits<T, U>::Promote
268 typename PromoteTraits<T, U>::Promote p(coeffs_[order_]);
269 for(
int i = order_ - 1; i >= 0; --i)
271 p = v * p + coeffs_[i];
301 for(
unsigned int i = 1; i <= order_; ++i)
303 coeffs_[i-1] = double(i)*coeffs_[i];
314 vigra_precondition(order_ > 0,
315 "PolynomialView<T>::deflate(): cannot deflate 0th order polynomial.");
326 forwardBackwardDeflate(v);
329 deflate(v, multiplicity-1);
336 for(
int i = order_-1; i > 0; --i)
338 coeffs_[i] += v * coeffs_[i+1];
348 unsigned int order2 = order_ / 2;
349 T tmp = coeffs_[order_];
350 for(
unsigned int i = order_-1; i >= order2; --i)
354 tmp = tmp1 + v * tmp;
358 for(
unsigned int i = 1; i < order2; ++i)
360 coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
371 for(
unsigned int i = 1; i < order_; ++i)
373 coeffs_[i] = v * (coeffs_[i] - coeffs_[i-1]);
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)
390 coeffs_[i] += a * coeffs_[i+1] + b*coeffs_[i+2];
400 while(
std::abs(coeffs_[order_]) <= epsilon && order_ > 0)
408 for(
unsigned int i = 0; i<order_; ++i)
409 coeffs_[i] /= coeffs_[order_];
410 coeffs_[order_] = T(1.0);
417 Real p0 =
abs(coeffs_[0]), po =
abs(coeffs_[order_]);
418 Real
norm = (p0 > 0.0)
421 for(
unsigned int i = 0; i<=order_; ++i)
465 polynomial_(
order + 1, T())
467 this->setCoeffs(&polynomial_[0],
order);
476 this->setCoeffs(&polynomial_[0], p.
order());
481 template <
class ITER>
484 polynomial_(i, i + order + 1)
486 this->setCoeffs(&polynomial_[0], order);
494 template <
class ITER>
497 polynomial_(i, i + order + 1)
499 this->setCoeffs(&polynomial_[0], order);
509 polynomial_.swap(tmp);
510 this->setCoeffs(&polynomial_[0], p.
order());
511 this->epsilon_ = p.epsilon_;
573 template <
unsigned int MAXORDER,
class T>
599 vigra_precondition(
order <= MAXORDER,
600 "StaticPolynomial(): order exceeds MAXORDER.");
601 std::fill_n(polynomial_,
order+1, T());
602 this->setCoeffs(polynomial_,
order);
610 std::copy(p.
begin(), p.
end(), polynomial_);
611 this->setCoeffs(polynomial_, p.
order());
617 template <
class ITER>
621 vigra_precondition(order <= MAXORDER,
622 "StaticPolynomial(): order exceeds MAXORDER.");
623 std::copy(i, i + order + 1, polynomial_);
624 this->setCoeffs(polynomial_, order);
632 template <
class ITER>
636 vigra_precondition(order <= MAXORDER,
637 "StaticPolynomial(): order exceeds MAXORDER.");
638 std::copy(i, i + order + 1, polynomial_);
639 this->setCoeffs(polynomial_, order);
648 std::copy(p.
begin(), p.
end(), polynomial_);
649 this->setCoeffs(polynomial_, p.
order());
650 this->epsilon_ = p.epsilon_;
688 void setOrder(
unsigned int order)
690 vigra_precondition(order <= MAXORDER,
691 "taticPolynomial::setOrder(): order exceeds MAXORDER.");
692 this->order_ =
order;
696 T polynomial_[MAXORDER+1];
706 std::complex<T> complexDiv(std::complex<T>
const & a, std::complex<T>
const & b)
708 const double abs_breal = b.real() < 0 ? -b.real() : b.real();
709 const double abs_bimag = b.imag() < 0 ? -b.imag() : b.imag();
711 if (abs_breal >= abs_bimag)
714 if (abs_breal == 0.0)
716 return std::complex<T>(a.real() / abs_breal, a.imag() / abs_breal);
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);
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);
737 std::complex<T> deleteBelowEpsilon(std::complex<T>
const & x,
double eps)
740 ? std::complex<T>(x.real())
742 ? std::complex<T>(NumericTraits<T>::zero(), x.imag())
746 template <
class POLYNOMIAL>
747 typename POLYNOMIAL::value_type
748 laguerreStartingGuess(POLYNOMIAL
const & p)
750 double N = p.order();
751 typename POLYNOMIAL::value_type centroid = -p[p.order()-1] / N / p[p.order()];
753 return centroid + dist;
756 template <
class POLYNOMIAL,
class Complex>
757 int laguerre1Root(POLYNOMIAL
const & p, Complex & x,
unsigned int multiplicity)
759 typedef typename NumericTraits<Complex>::ValueType Real;
761 static double frac[] = {0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
764 double N = p.order();
765 double eps = p.epsilon(),
768 if(multiplicity == 0)
769 x = laguerreStartingGuess(p);
771 bool mayTryDerivative =
true;
773 for(count = 0; count < maxiter; ++count)
777 Complex p0(p[p.order()]);
782 for(
int i = p.order()-1; i >= 0; --i)
796 Complex g = complexDiv(p1, p0);
798 Complex h = g2 - complexDiv(p2, p0);
803 (
std::abs(N * complexDiv(h, g2) - 1.0) + 1.0) + 0.5);
810 if(mayTryDerivative && multiplicity > 1 && ap0 < eps2)
813 int derivativeMultiplicity = laguerre1Root(p.getDerivative(), x1, multiplicity-1);
818 return derivativeMultiplicity + 1;
823 mayTryDerivative =
false;
834 dx = complexDiv(Complex(N) , gp);
851 x = x - frac[(count+1)/10] * dx;
853 return count < maxiter ?
858 template <
class Real>
859 struct PolynomialRootCompare
863 PolynomialRootCompare(Real eps)
868 bool operator()(T
const & l, T
const & r)
871 ? l.imag() < r.imag()
872 : l.real() < r.real();
938 template <
class POLYNOMIAL,
class VECTOR>
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;
946 double eps = poriginal.epsilon();
948 WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
953 Complex x = detail::laguerreStartingGuess(p);
955 unsigned int multiplicity = 1;
956 bool triedConjugate =
false;
965 multiplicity = detail::laguerre1Root(p, x, multiplicity);
967 if(multiplicity == 0)
971 if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
973 x = detail::deleteBelowEpsilon(x, eps);
981 triedConjugate =
false;
986 if(x.imag() != 0.0 && !triedConjugate)
990 triedConjugate =
true;
996 triedConjugate =
false;
997 x = detail::laguerreStartingGuess(p);
1010 if((
conj(b)*b2).real() >= 0.0)
1011 q = -0.5 * (b + b2);
1013 q = -0.5 * (b - b2);
1014 x = detail::complexDiv(q, a);
1016 detail::laguerre1Root(poriginal, x, 1);
1017 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1018 x = detail::complexDiv(c, q);
1020 detail::laguerre1Root(poriginal, x, 1);
1021 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1023 else if(p.order() == 1)
1025 x = detail::complexDiv(-p[0], p[1]);
1027 detail::laguerre1Root(poriginal, x, 1);
1028 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1030 std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
1034 template <
class POLYNOMIAL,
class VECTOR>
1076 template <
class POLYNOMIAL,
class VECTOR>
1079 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1083 for(
unsigned int i = 0; i < croots.
size(); ++i)
1084 if(croots[i].imag() == 0.0)
1085 roots.push_back(croots[i].real());
1089 template <
class POLYNOMIAL,
class VECTOR>
1103 ostream & operator<<(ostream & o, vigra::PolynomialView<T>
const & p)
1105 for(
unsigned int k=0; k < p.order(); ++k)
1113 #endif // VIGRA_POLYNOMIAL_HXX