00001 #ifndef vnl_rnpoly_solve_h_ 00002 #define vnl_rnpoly_solve_h_ 00003 00004 //: 00005 // \file 00006 // \brief Solves for roots of system of real polynomials 00007 // \author Marc Pollefeys, ESAT-VISICS, K.U.Leuven, 12-08-97 00008 // 00009 // Modifications 00010 // Peter Vanroose, 20 Oct 1999: implementation simplified through "cmplx" 00011 // class for doing complex arithmetic. 00012 // dac (Manchester) 28/03/2001: tidied up documentation 00013 // 00014 00015 #include <vnl/vnl_vector.h> 00016 #include <vnl/vnl_real_npolynomial.h> 00017 #include <vcl_vector.h> 00018 00019 //: Solves for roots of system of real polynomials 00020 // Calculates all the roots of a system of N polynomials in N variables 00021 // through continuation. 00022 // Adapted from the PARALLEL CONTINUATION algorithm , written by Darrell 00023 // Stam, 1991, and further improved by Kriegman and Ponce, 1992. 00024 // 00025 00026 00027 #ifdef static 00028 # error "grr!!" 00029 #endif 00030 00031 //: Solves for roots of system of real polynomials 00032 00033 class vnl_rnpoly_solve { 00034 public: 00035 #ifndef _WIN32 00036 static const unsigned int M = 11; // Maximum dimension of problem 00037 static const unsigned int T = 2500; // Max. number of terms in a polynomial 00038 #else 00039 enum { M = 11 }; // Maximum dimension of problem 00040 enum { T = 2500 }; // Maximum number of terms in a polynomial 00041 #endif 00042 00043 // Constructor--------------------------------------------------------------- 00044 00045 //: The constructor already does all the calculations 00046 inline vnl_rnpoly_solve(vcl_vector<vnl_real_npolynomial*> 00047 const& ps) : ps_(ps) {compute();} 00048 ~vnl_rnpoly_solve(); 00049 00050 // Operations---------------------------------------------------------------- 00051 00052 //: Array of real parts of roots 00053 inline vcl_vector<vnl_vector<double>*> real() { return r_; } 00054 00055 //: Array of imaginary parts of roots 00056 inline vcl_vector<vnl_vector<double>*> imag() { return i_; } 00057 00058 //: Return real roots only. Roots are real if the absolute value 00059 // of their imaginary part is less than the optional argument tol, 00060 // which defaults to 1e-12 [untested] 00061 vcl_vector<vnl_vector<double>*> realroots(double tol = 1e-12); 00062 00063 // Computations-------------------------------------------------------------- 00064 00065 private: 00066 //: Compute roots using continuation algorithm. 00067 bool compute(); 00068 00069 int Read_Input(int ideg[M], int terms[M], 00070 int polyn[M][T][M], double coeff[M][T]); 00071 00072 // Data Members-------------------------------------------------------------- 00073 vcl_vector<vnl_real_npolynomial*> ps_; // the input 00074 vcl_vector<vnl_vector<double>*> r_; // the output (real part) 00075 vcl_vector<vnl_vector<double>*> i_; // the output (imaginary part) 00076 }; 00077 00078 #endif // vnl_rnpoly_solve_h_