poly.cpp
Go to the documentation of this file.
#include <itpp/base/itcompat.h>
#include <itpp/signal/poly.h>
#include <itpp/base/converters.h>
#include <itpp/base/algebra/eigen.h>
#include <itpp/base/specmat.h>
#include <itpp/base/matfunc.h>
namespace itpp
{
void poly(const vec &r, vec &p)
{
int n = r.size();
p.set_size(n + 1, false);
p.zeros();
p(0) = 1.0;
for (int i = 0; i < n; i++)
p.set_subvector(1, p(1, i + 1) - r(i)*p(0, i));
}
void poly(const cvec &r, cvec &p)
{
int n = r.size();
p.set_size(n + 1, false);
p.zeros();
p(0) = 1.0;
for (int i = 0; i < n; i++)
p.set_subvector(1, p(1, i + 1) - r(i)*p(0, i));
}
void roots(const vec &p, cvec &r)
{
int n = p.size(), m, l;
ivec f = find(p != 0.0);
m = f.size();
vec v = p;
mat A;
if (m > 0 && n > 1) {
v = v(f(0), f(m - 1));
l = v.size();
if (l > 1) {
A = diag(ones(l - 2), -1);
A.set_row(0, -v(1, l - 1) / v(0));
r = eig(A);
cvec d;
cmat V;
eig(A, d , V);
if (f(m - 1) < n)
r = concat(r, zeros_c(n - f(m - 1) - 1));
}
else {
r.set_size(n - f(m - 1) - 1, false);
r.zeros();
}
}
else
r.set_size(0, false);
}
void roots(const cvec &p, cvec &r)
{
int n = p.size(), m, l;
ivec f;
// find all non-zero elements
for (int i = 0; i < n; i++)
if (p(i) != 0.0)
f = concat(f, i);
m = f.size();
cvec v = p;
cmat A;
if (m > 0 && n > 1) {
v = v(f(0), f(m - 1));
l = v.size();
if (l > 1) {
A = diag(ones_c(l - 2), -1);
A.set_row(0, -v(1, l - 1) / v(0));
r = eig(A);
if (f(m - 1) < n)
r = concat(r, zeros_c(n - f(m - 1) - 1));
}
else {
r.set_size(n - f(m - 1) - 1, false);
r.zeros();
}
}
else
r.set_size(0, false);
}
vec polyval(const vec &p, const vec &x)
{
it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
vec out(x.size());
out = p(0);
for (int i = 1; i < p.size(); i++)
out = p(i) + elem_mult(x, out);
return out;
}
cvec polyval(const vec &p, const cvec &x)
{
it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
cvec out(x.size());
out = p(0);
for (int i = 1; i < p.size(); i++)
out = std::complex<double>(p(i)) + elem_mult(x, out);
return out;
}
cvec polyval(const cvec &p, const vec &x)
{
it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
cvec out(x.size());
out = p(0);
for (int i = 1; i < p.size(); i++)
out = std::complex<double>(p(i)) + elem_mult(to_cvec(x), out);
return out;
}
cvec polyval(const cvec &p, const cvec &x)
{
it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
cvec out(x.size());
out = p(0);
for (int i = 1; i < p.size(); i++)
out = p(i) + elem_mult(x, out);
return out;
}
double cheb(int n, double x)
{
it_assert((n >= 0), "cheb(): need a non-negative order n!");
if (x < 1.0 && x > -1.0) {
return std::cos(n * std::acos(x));
}
else if (x <= -1) {
return (is_even(n) ? std::cosh(n * ::acosh(-x))
: -std::cosh(n * ::acosh(-x)));
}
return std::cosh(n * ::acosh(x));
}
vec cheb(int n, const vec &x)
{
it_assert_debug(x.size() > 0, "cheb(): empty vector");
vec out(x.size());
for (int i = 0; i < x.size(); ++i) {
out(i) = cheb(n, x(i));
}
return out;
}
mat cheb(int n, const mat &x)
{
it_assert_debug((x.rows() > 0) && (x.cols() > 0), "cheb(): empty matrix");
mat out(x.rows(), x.cols());
for (int i = 0; i < x.rows(); ++i) {
for (int j = 0; j < x.cols(); ++j) {
out(i, j) = cheb(n, x(i, j));
}
}
return out;
}
} // namespace itpp
itpp::concat
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
itpp::diag
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
eigen.h
Definitions of eigenvalue decomposition functions.
itpp::ones_c
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
itpp::is_even
bool is_even(int x)
Return true if x is an even integer.
Definition: misc.h:122
it_error_if
#define it_error_if(t, s)
Abort if t is true.
Definition: itassert.h:117
itpp::cosh
vec cosh(const vec &x)
Cosine hyperbolic function.
Definition: trig_hyp.h:93
itcompat.h
IT++ compatibility types and functions.
itpp::find
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition: specmat.cpp:40
itpp
itpp namespace
Definition: itmex.h:36
itpp::elem_mult
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
it_assert_debug
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
itpp::eig
bool eig(const mat &A, cvec &d, cmat &V)
Calculates the eigenvalues and eigenvectors of a real non-symmetric matrix.
Definition: eigen.cpp:277
matfunc.h
Various functions on vectors and matrices - header file.
itpp::acosh
vec acosh(const vec &x)
Inverse cosine hyperbolic function.
Definition: trig_hyp.cpp:40
itpp::to_cvec
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition: converters.h:107
specmat.h
Definitions of special vectors and matrices.
itpp::cheb
double cheb(int n, double x)
Chebyshev polynomial of the first kind.
Definition: poly.cpp:195
poly.h
Polynomial functions.
itpp::zeros_c
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
itpp::roots
void roots(const vec &p, cvec &r)
Calculate the roots of the polynomial.
Definition: poly.cpp:66
itpp::cos
vec cos(const vec &x)
Cosine function.
Definition: trig_hyp.h:58
itpp::polyval
vec polyval(const vec &p, const vec &x)
Evaluate polynomial.
Definition: poly.cpp:135
converters.h
Definitions of converters between different vector and matrix types.
itpp::poly
void poly(const vec &r, vec &p)
Create a polynomial of the given roots.
Definition: poly.cpp:40
itpp::ones
ITPP_EXPORT vec ones(int size)
A float vector of ones.
itpp::acos
vec acos(const vec &x)
Inverse cosine function.
Definition: trig_hyp.h:70
it_assert
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94

Generated on Mon Mar 23 2020 06:18:03 for IT++ by Doxygen 1.8.17