1 // This is core/vnl/algo/vnl_rpoly_roots.h
2 #ifndef vnl_rpoly_roots_h_
3 #define vnl_rpoly_roots_h_
4 //:
5 //  \file
6 //  \brief Finds roots of a real polynomial
7 //  \author  Andrew W. Fitzgibbon, Oxford RRG
8 //  \date   06 Aug 96
9 //
10 // \verbatim
11 //  Modifications
12 //  23 may 97, Peter Vanroose - "NO_COMPLEX" option added (until "complex" type is standardised)
13 //  dac (Manchester) 28/03/2001: tidied up documentation
14 //  Joris Van den Wyngaerd - June 2001 - impl for vnl_real_polynomial constr added
15 //   Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line
16 //  \endverbatim
17 
18 #include <complex>
19 #ifdef _MSC_VER
20 #  include <vcl_msvc_warnings.h>
21 #endif
22 #include <vnl/vnl_vector.h>
23 #include <vnl/algo/vnl_algo_export.h>
24 
25 class vnl_real_polynomial;
26 
27 //: Find the roots of a real polynomial.
28 //  Uses algorithm 493 from
29 //  ACM Trans. Math. Software - the Jenkins-Traub algorithm, described
30 //  by Numerical Recipes under "Other sure-fire techniques" as
31 //  "practically a standard in black-box polynomial rootfinders".
32 //  (See M.A. Jenkins, ACM TOMS 1 (1975) pp. 178-189.).
33 //
34 //  This class is not very const-correct as it is intended as a compute object
35 //  rather than a data object.
36 
37 class VNL_ALGO_EXPORT vnl_rpoly_roots
38 {
39  public:
40 // Constructors/Destructors--------------------------------------------------
41 
42   //: The constructor calculates the roots.
43   // This is the most efficient interface
44   // as all the result variables are initialized to the correct size.
45   // The polynomial is $ a[0] x^d + a[1] x^{d-1} + \cdots + a[d] = 0 $.
46   //
47   // Note that if the routine fails, not all roots will be found.  In this case,
48   // the "realroots" and "roots" functions will return fewer than n roots.
49 
50   vnl_rpoly_roots(const vnl_vector<double>& a);
51 
52   //: Calculate roots of a vnl_real_polynomial. Same comments apply.
53   vnl_rpoly_roots(const vnl_real_polynomial& poly);
54 
55   // Operations----------------------------------------------------------------
56 
57   //: Return i'th complex root
58   std::complex<double> operator [] (int i) const { return {r_[i], i_[i]}; }
59 
60   //: Complex vector of all roots.
61   vnl_vector<std::complex<double> > roots() const;
62 
63   //: Real part of root I.
real(int i)64   const double& real(int i) const { return r_[i]; }
65 
66   //: Imaginary part of root I.
imag(int i)67   const double& imag(int i) const { return i_[i]; }
68 
69   //: Vector of real parts of roots
real()70   vnl_vector<double>& real() { return r_; }
71 
72   //: Vector of imaginary parts of roots
imag()73   vnl_vector<double>& imag() { return i_; }
74 
75   //: Return real roots only.
76   //  Roots are real if the absolute value of their imaginary part is less than
77   //  the optional argument TOL. TOL defaults to 1e-12 [untested]
78   vnl_vector<double> realroots(double tol = 1e-12) const;
79 
80   // Computations--------------------------------------------------------------
81 
82   //: Compute roots using Jenkins-Traub algorithm.
83   bool compute();
84 
85   //: Compute roots using QR decomposition of companion matrix. [unimplemented]
86   bool compute_qr();
87 
88   //: Compute roots using Laguerre algorithm. [unimplemented]
89   bool compute_laguerre();
90 
91  protected:
92   // Data Members--------------------------------------------------------------
93   vnl_vector<double> coeffs_;
94 
95   vnl_vector<double> r_;
96   vnl_vector<double> i_;
97 
98   int num_roots_found_;
99 };
100 
101 #endif // vnl_rpoly_roots_h_
102