37 #ifndef VIGRA_POLYNOMIAL_HXX
38 #define VIGRA_POLYNOMIAL_HXX
45 #include "mathutil.hxx"
46 #include "numerictraits.hxx"
47 #include "array_vector.hxx"
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)
444 :
public PolynomialView<T>
446 typedef PolynomialView<T> BaseType;
450 typedef Polynomial<Real> RealPolynomial;
451 typedef Polynomial<Complex> ComplexPolynomial;
453 typedef T value_type;
454 typedef T * iterator;
455 typedef T
const * const_iterator;
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_;
521 res.differentiate(n);
573 template <
unsigned int MAXORDER,
class T>
574 class StaticPolynomial
575 :
public PolynomialView<T>
577 typedef PolynomialView<T> BaseType;
582 typedef StaticPolynomial<MAXORDER, Real> RealPolynomial;
583 typedef StaticPolynomial<MAXORDER, Complex> ComplexPolynomial;
585 typedef T value_type;
586 typedef T * iterator;
587 typedef T
const * const_iterator;
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_;
660 res.differentiate(n);
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()];
752 double dist = VIGRA_CSTD::pow(
std::abs(p(centroid) / p[p.order()]), 1.0 / N);
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 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::Real Real;
942 typedef typename POLYNOMIAL::Complex Complex;
943 typedef typename POLYNOMIAL::ComplexPolynomial WorkPolynomial;
945 double eps = poriginal.epsilon();
947 WorkPolynomial p(poriginal.begin(), poriginal.order(), eps);
952 Complex x = detail::laguerreStartingGuess(p);
954 unsigned int multiplicity = 1;
955 bool triedConjugate =
false;
964 multiplicity = detail::laguerre1Root(p, x, multiplicity);
966 if(multiplicity == 0)
970 if(polishRoots && !detail::laguerre1Root(poriginal, x, multiplicity))
972 x = detail::deleteBelowEpsilon(x, eps);
980 triedConjugate =
false;
985 if(x.imag() != 0.0 && !triedConjugate)
989 triedConjugate =
true;
995 triedConjugate =
false;
996 x = detail::laguerreStartingGuess(p);
1010 q = -0.5 * (b + b2);
1012 q = -0.5 * (b - b2);
1013 x = detail::complexDiv(q, a);
1015 detail::laguerre1Root(poriginal, x, 1);
1016 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1017 x = detail::complexDiv(c, q);
1019 detail::laguerre1Root(poriginal, x, 1);
1020 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1022 else if(p.order() == 1)
1024 x = detail::complexDiv(-p[0], p[1]);
1026 detail::laguerre1Root(poriginal, x, 1);
1027 roots.push_back(detail::deleteBelowEpsilon(x, eps));
1029 std::sort(roots.begin(), roots.end(), detail::PolynomialRootCompare<Real>(eps));
1033 template <
class POLYNOMIAL,
class VECTOR>
1075 template <
class POLYNOMIAL,
class VECTOR>
1078 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
1082 for(
unsigned int i = 0; i < croots.
size(); ++i)
1083 if(croots[i].
imag() == 0.0)
1084 roots.push_back(croots[i].
real());
1088 template <
class POLYNOMIAL,
class VECTOR>
1102 ostream & operator<<(ostream & o, vigra::PolynomialView<T>
const & p)
1104 for(
unsigned int k=0; k < p.order(); ++k)
1112 #endif // VIGRA_POLYNOMIAL_HXX
Polynomial< T > getDeflated(Real r) const
Definition: polynomial.hxx:529
StaticPolynomial(ITER i, unsigned int order)
Definition: polynomial.hxx:618
StaticPolynomial & operator=(StaticPolynomial const &p)
Definition: polynomial.hxx:644
Definition: polynomial.hxx:51
const_iterator end() const
Definition: polynomial.hxx:220
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
R imag(const FFTWComplex< R > &a)
imaginary part
Definition: fftw3.hxx:1023
double epsilon() const
Definition: polynomial.hxx:236
Polynomial & operator=(Polynomial const &p)
Definition: polynomial.hxx:504
T value_type
Definition: polynomial.hxx:79
bool polynomialRoots(POLYNOMIAL const &poriginal, VECTOR &roots, bool polishRoots)
Definition: polynomial.hxx:939
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
Polynomial(Polynomial const &p)
Definition: polynomial.hxx:472
void setEpsilon(double eps)
Definition: polynomial.hxx:242
Definition: polynomial.hxx:52
T const * const_iterator
Definition: polynomial.hxx:100
T * iterator
Definition: polynomial.hxx:96
R real(const FFTWComplex< R > &a)
real part
Definition: fftw3.hxx:1016
PolynomialView(T *coeffs, unsigned int order, double epsilon=1.0e-14)
Definition: polynomial.hxx:120
Polynomial< Complex > getDeflated(Complex const &r) const
Definition: polynomial.hxx:542
T & operator[](unsigned int i)
Access the coefficient of x^i.
Definition: polynomial.hxx:127
const_iterator begin() const
Definition: polynomial.hxx:215
NumericTraits< RealPromote >::ValueType Real
Definition: polynomial.hxx:88
StaticPolynomial(ITER i, unsigned int order, double epsilon)
Definition: polynomial.hxx:633
Polynomial< T > getDerivative(unsigned int n=1) const
Definition: polynomial.hxx:518
FixedPoint< 0, FracBits > frac(FixedPoint< IntBits, FracBits > v)
fractional part.
Definition: fixedpoint.hxx:650
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
void deflate(T const &r, unsigned int multiplicity=1)
Definition: polynomial.hxx:312
iterator end()
Definition: polynomial.hxx:210
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
StaticPolynomial(StaticPolynomial const &p)
Definition: polynomial.hxx:607
T const & operator[](unsigned int i) const
Access the coefficient of x^i.
Definition: polynomial.hxx:131
unsigned int order() const
Definition: polynomial.hxx:230
StaticPolynomial< MAXORDER, Complex > getDeflated(Complex const &r) const
Definition: polynomial.hxx:681
bool polynomialRealRoots(POLYNOMIAL const &p, VECTOR &roots, bool polishRoots)
Definition: polynomial.hxx:1076
Polynomial(ITER i, unsigned int order, double epsilon)
Definition: polynomial.hxx:495
StaticPolynomial(unsigned int order=0, double epsilon=1.0e-14)
Definition: polynomial.hxx:596
void forwardDeflate(T const &v)
Definition: polynomial.hxx:334
NumericTraits< RealPromote >::ComplexPromote Complex
Definition: polynomial.hxx:92
StaticPolynomial getDeflated(Real r) const
Definition: polynomial.hxx:668
unsigned int size() const
Definition: polynomial.hxx:225
iterator begin()
Definition: polynomial.hxx:205
NumericTraits< T >::RealPromote RealPromote
Definition: polynomial.hxx:84
Polynomial(ITER i, unsigned int order)
Definition: polynomial.hxx:482
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point equality.
Definition: mathutil.hxx:1638
StaticPolynomial getDerivative(unsigned int n=1) const
Definition: polynomial.hxx:657
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
void differentiate(unsigned int n=1)
Definition: polynomial.hxx:292
void backwardDeflate(T v)
Definition: polynomial.hxx:367
void minimizeOrder(double epsilon=0.0)
Definition: polynomial.hxx:398
Definition: polynomial.hxx:73
size_type size() const
Definition: array_vector.hxx:358
void normalize()
Definition: polynomial.hxx:406
PromoteTraits< T, V >::Promote operator()(V v) const
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
Polynomial(unsigned int order=0, double epsilon=1.0e-14)
Definition: polynomial.hxx:463
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
void forwardBackwardDeflate(T v)
Definition: polynomial.hxx:346
void deflateConjugatePair(Complex const &v)
Definition: polynomial.hxx:380
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616