00001 00029 #include <itpp/base/math/elem_math.h> 00030 #include <itpp/base/itcompat.h> 00031 00032 00033 namespace itpp 00034 { 00035 00036 vec sqr(const cvec &data) 00037 { 00038 vec temp(data.length()); 00039 for (int i = 0; i < data.length(); i++) 00040 temp(i) = sqr(data(i)); 00041 return temp; 00042 } 00043 00044 mat sqr(const cmat &data) 00045 { 00046 mat temp(data.rows(), data.cols()); 00047 for (int i = 0; i < temp.rows(); i++) { 00048 for (int j = 0; j < temp.cols(); j++) { 00049 temp(i, j) = sqr(data(i, j)); 00050 } 00051 } 00052 return temp; 00053 } 00054 00055 vec abs(const cvec &data) 00056 { 00057 vec temp(data.length()); 00058 00059 for (int i = 0;i < data.length();i++) 00060 temp[i] = std::abs(data[i]); 00061 00062 return temp; 00063 } 00064 00065 mat abs(const cmat &data) 00066 { 00067 mat temp(data.rows(), data.cols()); 00068 00069 for (int i = 0;i < temp.rows();i++) { 00070 for (int j = 0;j < temp.cols();j++) { 00071 temp(i, j) = std::abs(data(i, j)); 00072 } 00073 } 00074 00075 return temp; 00076 } 00077 00078 // Deprecated gamma function. Will be changed to tgamma(). 00079 double gamma(double x) { return tgamma(x); } 00080 vec gamma(const vec &x) { return apply_function<double>(tgamma, x); } 00081 mat gamma(const mat &x) { return apply_function<double>(tgamma, x); } 00082 00083 // Calculates factorial coefficient for index <= 170. 00084 double fact(int index) 00085 { 00086 it_error_if(index > 170, "fact(int index): Function overflows if index > 170."); 00087 it_error_if(index < 0, "fact(int index): index must be non-negative integer"); 00088 double prod = 1; 00089 for (int i = 1; i <= index; i++) 00090 prod *= static_cast<double>(i); 00091 return prod; 00092 } 00093 00094 // Calculates binomial coefficient "n over k". 00095 double binom(int n, int k) 00096 { 00097 it_assert(k <= n, "binom(n, k): k can not be larger than n"); 00098 it_assert((n >= 0) && (k >= 0), "binom(n, k): n and k must be non-negative integers"); 00099 k = ((n - k) < k) ? n - k : k; 00100 00101 double out = 1.0; 00102 for (int i = 1; i <= k; ++i) { 00103 out *= (i + n - k); 00104 out /= i; 00105 } 00106 return out; 00107 } 00108 00109 // Calculates binomial coefficient "n over k". 00110 int binom_i(int n, int k) 00111 { 00112 it_assert(k <= n, "binom_i(n, k): k can not be larger than n"); 00113 it_assert((n >= 0) && (k >= 0), "binom_i(n, k): n and k must be non-negative integers"); 00114 k = ((n - k) < k) ? n - k : k; 00115 00116 int out = 1; 00117 for (int i = 1; i <= k; ++i) { 00118 out *= (i + n - k); 00119 out /= i; 00120 } 00121 return out; 00122 } 00123 00124 // Calculates the base 10-logarithm of the binomial coefficient "n over k". 00125 double log_binom(int n, int k) 00126 { 00127 it_assert(k <= n, "log_binom(n, k): k can not be larger than n"); 00128 it_assert((n >= 0) && (k >= 0), "log_binom(n, k): n and k must be non-negative integers"); 00129 k = ((n - k) < k) ? n - k : k; 00130 00131 double out = 0.0; 00132 for (int i = 1; i <= k; i++) 00133 out += log10(static_cast<double>(i + n - k)) 00134 - log10(static_cast<double>(i)); 00135 00136 return out; 00137 } 00138 00139 // Calculates the greatest common divisor 00140 int gcd(int a, int b) 00141 { 00142 it_assert((a >= 0) && (b >= 0), "gcd(a, b): a and b must be non-negative integers"); 00143 int v, u, t, q; 00144 00145 u = a; 00146 v = b; 00147 while (v > 0) { 00148 q = u / v; 00149 t = u - v * q; 00150 u = v; 00151 v = t; 00152 } 00153 return u; 00154 } 00155 00156 00157 vec real(const cvec &data) 00158 { 00159 vec temp(data.length()); 00160 00161 for (int i = 0;i < data.length();i++) 00162 temp[i] = data[i].real(); 00163 00164 return temp; 00165 } 00166 00167 mat real(const cmat &data) 00168 { 00169 mat temp(data.rows(), data.cols()); 00170 00171 for (int i = 0;i < temp.rows();i++) { 00172 for (int j = 0;j < temp.cols();j++) { 00173 temp(i, j) = data(i, j).real(); 00174 } 00175 } 00176 00177 return temp; 00178 } 00179 00180 vec imag(const cvec &data) 00181 { 00182 vec temp(data.length()); 00183 00184 for (int i = 0;i < data.length();i++) 00185 temp[i] = data[i].imag(); 00186 return temp; 00187 } 00188 00189 mat imag(const cmat &data) 00190 { 00191 mat temp(data.rows(), data.cols()); 00192 00193 for (int i = 0;i < temp.rows();i++) { 00194 for (int j = 0;j < temp.cols();j++) { 00195 temp(i, j) = data(i, j).imag(); 00196 } 00197 } 00198 00199 return temp; 00200 } 00201 00202 vec arg(const cvec &data) 00203 { 00204 vec temp(data.length()); 00205 00206 for (int i = 0;i < data.length();i++) 00207 temp[i] = std::arg(data[i]); 00208 00209 return temp; 00210 } 00211 00212 mat arg(const cmat &data) 00213 { 00214 mat temp(data.rows(), data.cols()); 00215 00216 for (int i = 0;i < temp.rows();i++) { 00217 for (int j = 0;j < temp.cols();j++) { 00218 temp(i, j) = std::arg(data(i, j)); 00219 } 00220 } 00221 00222 return temp; 00223 } 00224 00225 #ifdef _MSC_VER 00226 cvec conj(const cvec &x) 00227 { 00228 cvec temp(x.size()); 00229 00230 for (int i = 0; i < x.size(); i++) { 00231 temp(i) = std::conj(x(i)); 00232 } 00233 00234 return temp; 00235 } 00236 00237 cmat conj(const cmat &x) 00238 { 00239 cmat temp(x.rows(), x.cols()); 00240 00241 for (int i = 0; i < x.rows(); i++) { 00242 for (int j = 0; j < x.cols(); j++) { 00243 temp(i, j) = std::conj(x(i, j)); 00244 } 00245 } 00246 00247 return temp; 00248 } 00249 #endif 00250 00251 } // namespace itpp
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4