1 /*
2 fsm
3 */
4 #include <cassert>
5 #include "vnl_cpoly_roots.h"
6 #include <vnl/algo/vnl_complex_eigensystem.h>
7
compute(vnl_vector<std::complex<double>> const & a)8 void vnl_cpoly_roots::compute(vnl_vector<std::complex<double> > const &a)
9 {
10 // construct companion matrix
11 vnl_matrix<std::complex<double> > comp(N, N);
12 comp.fill(0);
13 for (unsigned i=0; i<N-1; ++i)
14 comp(i+1, i) = 1;
15 for (unsigned i=0; i<N; ++i)
16 comp(i, N-1) = -a[N-1-i];
17
18 // the eigenvalues of the companion matrix are the roots of the polynomial
19 solns = vnl_complex_eigensystem(comp,
20 false, // we only want
21 false).W; // the eigenvalues.
22 #ifdef DEBUG
23 std::cerr << "s = " << solns << '\n';
24 #endif
25 }
26
vnl_cpoly_roots(vnl_vector<std::complex<double>> const & a)27 vnl_cpoly_roots::vnl_cpoly_roots(vnl_vector<std::complex<double> > const & a)
28 : solns(a.size())
29 , N(a.size()) // degree
30 {
31 compute(a);
32 }
33
vnl_cpoly_roots(vnl_vector<double> const & a_real,vnl_vector<double> const & a_imag)34 vnl_cpoly_roots::vnl_cpoly_roots(vnl_vector<double> const & a_real,
35 vnl_vector<double> const & a_imag)
36 : solns(a_real.size())
37 , N(a_real.size()) // degree
38 {
39 assert(a_real.size() == a_imag.size());
40 vnl_vector<std::complex<double> > a(N);
41 for (unsigned i=0; i<N; ++i)
42 a[i] = std::complex<double>(a_real[i], a_imag[i]);
43
44 #ifdef DEBUG
45 std::cerr << "a = " << a << '\n';
46 #endif
47 compute(a);
48 }
49