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