IT++ Logo
poly.cpp
Go to the documentation of this file.
1 
29 #include <itpp/base/itcompat.h>
30 #include <itpp/signal/poly.h>
31 #include <itpp/base/converters.h>
33 #include <itpp/base/specmat.h>
34 #include <itpp/base/matfunc.h>
35 
36 
37 namespace itpp
38 {
39 
40 void poly(const vec &r, vec &p)
41 {
42  int n = r.size();
43 
44  p.set_size(n + 1, false);
45  p.zeros();
46  p(0) = 1.0;
47 
48  for (int i = 0; i < n; i++)
49  p.set_subvector(1, p(1, i + 1) - r(i)*p(0, i));
50 }
51 
52 void poly(const cvec &r, cvec &p)
53 {
54  int n = r.size();
55 
56  p.set_size(n + 1, false);
57  p.zeros();
58  p(0) = 1.0;
59 
60  for (int i = 0; i < n; i++)
61  p.set_subvector(1, p(1, i + 1) - r(i)*p(0, i));
62 }
63 
64 
65 
66 void roots(const vec &p, cvec &r)
67 {
68  int n = p.size(), m, l;
69  ivec f = find(p != 0.0);
70  m = f.size();
71  vec v = p;
72  mat A;
73 
74  if (m > 0 && n > 1) {
75  v = v(f(0), f(m - 1));
76  l = v.size();
77 
78  if (l > 1) {
79 
80  A = diag(ones(l - 2), -1);
81  A.set_row(0, -v(1, l - 1) / v(0));
82  r = eig(A);
83  cvec d;
84  cmat V;
85  eig(A, d , V);
86 
87  if (f(m - 1) < n)
88  r = concat(r, zeros_c(n - f(m - 1) - 1));
89  }
90  else {
91  r.set_size(n - f(m - 1) - 1, false);
92  r.zeros();
93  }
94  }
95  else
96  r.set_size(0, false);
97 }
98 
99 void roots(const cvec &p, cvec &r)
100 {
101  int n = p.size(), m, l;
102  ivec f;
103 
104  // find all non-zero elements
105  for (int i = 0; i < n; i++)
106  if (p(i) != 0.0)
107  f = concat(f, i);
108 
109 
110  m = f.size();
111  cvec v = p;
112  cmat A;
113 
114  if (m > 0 && n > 1) {
115  v = v(f(0), f(m - 1));
116  l = v.size();
117 
118  if (l > 1) {
119  A = diag(ones_c(l - 2), -1);
120  A.set_row(0, -v(1, l - 1) / v(0));
121  r = eig(A);
122  if (f(m - 1) < n)
123  r = concat(r, zeros_c(n - f(m - 1) - 1));
124  }
125  else {
126  r.set_size(n - f(m - 1) - 1, false);
127  r.zeros();
128  }
129  }
130  else
131  r.set_size(0, false);
132 }
133 
134 
135 vec polyval(const vec &p, const vec &x)
136 {
137  it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
138  it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
139 
140  vec out(x.size());
141 
142  out = p(0);
143 
144  for (int i = 1; i < p.size(); i++)
145  out = p(i) + elem_mult(x, out);
146 
147  return out;
148 }
149 
150 cvec polyval(const vec &p, const cvec &x)
151 {
152  it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
153  it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
154 
155  cvec out(x.size());
156 
157  out = p(0);
158 
159  for (int i = 1; i < p.size(); i++)
160  out = std::complex<double>(p(i)) + elem_mult(x, out);
161 
162  return out;
163 }
164 
165 cvec polyval(const cvec &p, const vec &x)
166 {
167  it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
168  it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
169 
170  cvec out(x.size());
171 
172  out = p(0);
173 
174  for (int i = 1; i < p.size(); i++)
175  out = std::complex<double>(p(i)) + elem_mult(to_cvec(x), out);
176 
177  return out;
178 }
179 
180 cvec polyval(const cvec &p, const cvec &x)
181 {
182  it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
183  it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
184 
185  cvec out(x.size());
186 
187  out = p(0);
188 
189  for (int i = 1; i < p.size(); i++)
190  out = p(i) + elem_mult(x, out);
191 
192  return out;
193 }
194 
195 double cheb(int n, double x)
196 {
197  it_assert((n >= 0), "cheb(): need a non-negative order n!");
198 
199  if (x < 1.0 && x > -1.0) {
200  return std::cos(n * std::acos(x));
201  }
202  else if (x <= -1) {
203  return (is_even(n) ? std::cosh(n * ::acosh(-x))
204  : -std::cosh(n * ::acosh(-x)));
205  }
206  return std::cosh(n * ::acosh(x));
207 }
208 
209 vec cheb(int n, const vec &x)
210 {
211  it_assert_debug(x.size() > 0, "cheb(): empty vector");
212 
213  vec out(x.size());
214  for (int i = 0; i < x.size(); ++i) {
215  out(i) = cheb(n, x(i));
216  }
217  return out;
218 }
219 
220 mat cheb(int n, const mat &x)
221 {
222  it_assert_debug((x.rows() > 0) && (x.cols() > 0), "cheb(): empty matrix");
223 
224  mat out(x.rows(), x.cols());
225  for (int i = 0; i < x.rows(); ++i) {
226  for (int j = 0; j < x.cols(); ++j) {
227  out(i, j) = cheb(n, x(i, j));
228  }
229  }
230  return out;
231 }
232 
233 } // namespace itpp
#define it_error_if(t, s)
Abort if t is true.
Definition: itassert.h:117
Various functions on vectors and matrices - header file.
vec cosh(const vec &x)
Cosine hyperbolic function.
Definition: trig_hyp.h:93
vec acos(const vec &x)
Inverse cosine function.
Definition: trig_hyp.h:70
Definitions of eigenvalue decomposition functions.
bool is_even(int x)
Return true if x is an even integer.
Definition: misc.h:122
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition: specmat.cpp:40
Mat< Num_T > elem_mult(const Mat< Num_T > &m1, const Mat< Num_T > &m2)
Element wise multiplication of two matrices.
Definition: mat.h:1582
Polynomial functions.
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
double cheb(int n, double x)
Chebyshev polynomial of the first kindChebyshev polynomials of the first kind can be defined as follo...
Definition: poly.cpp:195
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
Definitions of converters between different vector and matrix types.
vec acosh(const vec &x)
Inverse cosine hyperbolic function.
Definition: trig_hyp.cpp:40
vec polyval(const vec &p, const vec &x)
Evaluate polynomialEvaluate the polynomial p (of length at the points x The output is given by ...
Definition: poly.cpp:135
itpp namespace
Definition: itmex.h:36
ITPP_EXPORT vec ones(int size)
A float vector of ones.
IT++ compatibility types and functions.
bool eig(const mat &A, cvec &d, cmat &V)
Calculates the eigenvalues and eigenvectors of a real non-symmetric matrix.
Definition: eigen.cpp:277
Definitions of special vectors and matrices.
void roots(const vec &p, cvec &r)
Calculate the roots of the polynomialCalculate the roots r of the polynomial p.
Definition: poly.cpp:66
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition: converters.h:107
Mat< T > diag(const Vec< T > &v, const int K=0)
Create a diagonal matrix using vector v as its diagonal.
Definition: matfunc.h:557
vec cos(const vec &x)
Cosine function.
Definition: trig_hyp.h:58
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
void poly(const vec &r, vec &p)
Create a polynomial of the given rootsCreate a polynomial p with roots r.
Definition: poly.cpp:40
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition: array.h:486
SourceForge Logo

Generated on Sun Apr 10 2022 12:00:00 for IT++ by Doxygen 1.8.14