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
Generated on Wed Apr 18 11:45:33 2007 for IT++ by Doxygen 1.5.2