IT++ Logo Newcom Logo

poly.cpp

Go to the documentation of this file.
00001 
00033 #include <itpp/base/poly.h>
00034 #include <itpp/base/converters.h>
00035 #include <itpp/base/eigen.h>
00036 #include <itpp/base/specmat.h>
00037 #include <itpp/base/matfunc.h>
00038 
00039 
00040 namespace itpp {
00041 
00042   void poly(const vec &r, vec &p)
00043   {
00044     int n = r.size();
00045     
00046     p.set_size(n+1, false);
00047     p.zeros();
00048     p(0) = 1.0;
00049 
00050     for (int i=0; i<n; i++)
00051       p.set_subvector(1, i+1, p(1,i+1) - r(i)*p(0,i));
00052   }
00053 
00054   void poly(const cvec &r, cvec &p)
00055   {
00056     int n = r.size();
00057     
00058     p.set_size(n+1, false);
00059     p.zeros();
00060     p(0) = 1.0;
00061 
00062     for (int i=0; i<n; i++)
00063       p.set_subvector(1, i+1, p(1,i+1) - r(i)*p(0,i));
00064   }
00065 
00066 
00067 
00068   void roots(const vec &p, cvec &r)
00069   {
00070     int n = p.size(), m, l;
00071     ivec f = find(p != 0.0);
00072     m = f.size();
00073     vec v = p;
00074     mat A;
00075     
00076     if (m > 0 && n > 1) {
00077       v = v(f(0),f(m-1));
00078       l = v.size();
00079 
00080       if (l>1) {
00081         
00082         A = diag(ones(l-2), -1);
00083         A.set_row(0, -v(1,l-1)/v(0)); 
00084         r = eig(A);
00085         cvec d;
00086         cmat V;
00087         eig(A, d ,V);
00088 
00089         if (f(m-1) < n)
00090           r = concat(r, zeros_c(n-f(m-1)-1));
00091       } else {
00092         r.set_size(n-f(m-1)-1, false);
00093         r.zeros();
00094       }
00095     } else
00096       r.set_size(0, false);
00097   }
00098 
00099   void roots(const cvec &p, cvec &r)
00100   {
00101     int n = p.size(), m, l;
00102     ivec f;
00103 
00104     // find all non-zero elements
00105     for (int i=0; i<n; i++)
00106       if( p(i) != 0.0 )
00107         f = concat(f, i);
00108 
00109 
00110     m = f.size();
00111     cvec v = p;
00112     cmat A;
00113     
00114     if (m > 0 && n > 1) {
00115       v = v(f(0),f(m-1));
00116       l = v.size();
00117 
00118       if (l>1) {
00119         A = diag(ones_c(l-2), -1);
00120         A.set_row(0, -v(1,l-1)/v(0)); 
00121         r = eig(A);
00122         if (f(m-1) < n)
00123           r = concat(r, zeros_c(n-f(m-1)-1));
00124       } else {
00125         r.set_size(n-f(m-1)-1, false);
00126         r.zeros();
00127       }
00128     } else
00129       r.set_size(0, false);
00130   }
00131 
00132 
00133   vec polyval(const vec &p, const vec &x)
00134   {
00135     it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00136     it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00137 
00138     vec out(x.size());
00139 
00140     out = p(0);
00141 
00142     for (int i=1; i<p.size(); i++)
00143       out = p(i) + elem_mult(x, out);
00144 
00145     return out;
00146   }
00147 
00148   cvec polyval(const vec &p, const cvec &x)
00149   {
00150     it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00151     it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00152 
00153     cvec out(x.size());
00154 
00155     out = p(0);
00156 
00157     for (int i=1; i<p.size(); i++)
00158       out = std::complex<double>(p(i)) + elem_mult(x, out);
00159 
00160     return out;
00161   }
00162 
00163   cvec polyval(const cvec &p, const vec &x)
00164   {
00165     it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00166     it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00167 
00168     cvec out(x.size());
00169 
00170     out = p(0);
00171 
00172     for (int i=1; i<p.size(); i++)
00173       out = std::complex<double>(p(i)) + elem_mult(to_cvec(x), out);
00174 
00175     return out;
00176   }
00177 
00178   cvec polyval(const cvec &p, const cvec &x)
00179   {
00180     it_error_if(p.size() == 0, "polyval: size of polynomial is zero");
00181     it_error_if(x.size() == 0, "polyval: size of input value vector is zero");
00182 
00183     cvec out(x.size());
00184 
00185     out = p(0);
00186 
00187     for (int i=1; i<p.size(); i++)
00188       out = p(i) + elem_mult(x, out);
00189 
00190     return out;
00191   }
00192 
00193 
00194 
00195 } // namespace itpp
SourceForge Logo

Generated on Thu Aug 30 02:47:18 2007 for IT++ by Doxygen 1.5.3