00001 00030 #include <itpp/base/math/elem_math.h> 00031 00032 00033 #ifndef HAVE_TGAMMA 00034 // "True" gamma function 00035 double tgamma(double x) 00036 { 00037 double s = (2.50662827510730 + 190.9551718930764 / (x + 1) 00038 - 216.8366818437280 / (x + 2) + 60.19441764023333 00039 / (x + 3) - 3.08751323928546 / (x + 4) + 0.00302963870525 00040 / (x + 5) - 0.00001352385959072596 / (x + 6)) / x; 00041 if (s < 0) 00042 return -exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(-s)); 00043 else 00044 return exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(s)); 00045 } 00046 #endif 00047 00048 #if !defined(HAVE_LGAMMA) || (HAVE_DECL_SIGNGAM != 1) 00049 // The sign of the Gamma function is returned in the external integer 00050 // signgam declared in <math.h>. It is 1 when the Gamma function is positive 00051 // or zero, -1 when it is negative. However, MinGW definition of lgamma() 00052 // function does not use the global signgam variable. 00053 int signgam; 00054 // Logarithm of an absolute value of gamma function 00055 double lgamma(double x) 00056 { 00057 double gam = tgamma(x); 00058 signgam = (gam < 0) ? -1 : 1; 00059 return log(fabs(gam)); 00060 } 00061 #endif 00062 00063 #ifndef HAVE_CBRT 00064 // Cubic root 00065 double cbrt(double x) { return std::pow(x, 1.0 / 3.0); } 00066 #endif 00067 00068 00069 namespace itpp 00070 { 00071 00072 vec sqr(const cvec &data) 00073 { 00074 vec temp(data.length()); 00075 for (int i = 0; i < data.length(); i++) 00076 temp(i) = sqr(data(i)); 00077 return temp; 00078 } 00079 00080 mat sqr(const cmat &data) 00081 { 00082 mat temp(data.rows(), data.cols()); 00083 for (int i = 0; i < temp.rows(); i++) { 00084 for (int j = 0; j < temp.cols(); j++) { 00085 temp(i, j) = sqr(data(i, j)); 00086 } 00087 } 00088 return temp; 00089 } 00090 00091 vec abs(const cvec &data) 00092 { 00093 vec temp(data.length()); 00094 00095 for (int i = 0;i < data.length();i++) 00096 temp[i] = std::abs(data[i]); 00097 00098 return temp; 00099 } 00100 00101 mat abs(const cmat &data) 00102 { 00103 mat temp(data.rows(), data.cols()); 00104 00105 for (int i = 0;i < temp.rows();i++) { 00106 for (int j = 0;j < temp.cols();j++) { 00107 temp(i, j) = std::abs(data(i, j)); 00108 } 00109 } 00110 00111 return temp; 00112 } 00113 00114 // Calculates factorial coefficient for index <= 170. 00115 double fact(int index) 00116 { 00117 it_error_if(index > 170, "fact(int index): Function overflows if index > 170."); 00118 it_error_if(index < 0, "fact(int index): index must be non-negative integer"); 00119 double prod = 1; 00120 for (int i = 1; i <= index; i++) 00121 prod *= static_cast<double>(i); 00122 return prod; 00123 } 00124 00125 // Calculates binomial coefficient "n over k". 00126 double binom(int n, int k) 00127 { 00128 it_assert(k <= n, "binom(n, k): k can not be larger than n"); 00129 it_assert((n >= 0) && (k >= 0), "binom(n, k): n and k must be non-negative integers"); 00130 k = ((n - k) < k) ? n - k : k; 00131 00132 double out = 1.0; 00133 for (int i = 1; i <= k; ++i) { 00134 out *= (i + n - k); 00135 out /= i; 00136 } 00137 return out; 00138 } 00139 00140 // Calculates binomial coefficient "n over k". 00141 int binom_i(int n, int k) 00142 { 00143 it_assert(k <= n, "binom_i(n, k): k can not be larger than n"); 00144 it_assert((n >= 0) && (k >= 0), "binom_i(n, k): n and k must be non-negative integers"); 00145 k = ((n - k) < k) ? n - k : k; 00146 00147 int out = 1; 00148 for (int i = 1; i <= k; ++i) { 00149 out *= (i + n - k); 00150 out /= i; 00151 } 00152 return out; 00153 } 00154 00155 // Calculates the base 10-logarithm of the binomial coefficient "n over k". 00156 double log_binom(int n, int k) 00157 { 00158 it_assert(k <= n, "log_binom(n, k): k can not be larger than n"); 00159 it_assert((n >= 0) && (k >= 0), "log_binom(n, k): n and k must be non-negative integers"); 00160 k = ((n - k) < k) ? n - k : k; 00161 00162 double out = 0.0; 00163 for (int i = 1; i <= k; i++) 00164 out += log10(static_cast<double>(i + n - k)) 00165 - log10(static_cast<double>(i)); 00166 00167 return out; 00168 } 00169 00170 // Calculates the greatest common divisor 00171 int gcd(int a, int b) 00172 { 00173 it_assert((a >= 0) && (b >= 0), "gcd(a, b): a and b must be non-negative integers"); 00174 int v, u, t, q; 00175 00176 u = a; 00177 v = b; 00178 while (v > 0) { 00179 q = u / v; 00180 t = u - v * q; 00181 u = v; 00182 v = t; 00183 } 00184 return u; 00185 } 00186 00187 00188 vec real(const cvec &data) 00189 { 00190 vec temp(data.length()); 00191 00192 for (int i = 0;i < data.length();i++) 00193 temp[i] = data[i].real(); 00194 00195 return temp; 00196 } 00197 00198 mat real(const cmat &data) 00199 { 00200 mat temp(data.rows(), data.cols()); 00201 00202 for (int i = 0;i < temp.rows();i++) { 00203 for (int j = 0;j < temp.cols();j++) { 00204 temp(i, j) = data(i, j).real(); 00205 } 00206 } 00207 00208 return temp; 00209 } 00210 00211 vec imag(const cvec &data) 00212 { 00213 vec temp(data.length()); 00214 00215 for (int i = 0;i < data.length();i++) 00216 temp[i] = data[i].imag(); 00217 return temp; 00218 } 00219 00220 mat imag(const cmat &data) 00221 { 00222 mat temp(data.rows(), data.cols()); 00223 00224 for (int i = 0;i < temp.rows();i++) { 00225 for (int j = 0;j < temp.cols();j++) { 00226 temp(i, j) = data(i, j).imag(); 00227 } 00228 } 00229 00230 return temp; 00231 } 00232 00233 vec arg(const cvec &data) 00234 { 00235 vec temp(data.length()); 00236 00237 for (int i = 0;i < data.length();i++) 00238 temp[i] = std::arg(data[i]); 00239 00240 return temp; 00241 } 00242 00243 mat arg(const cmat &data) 00244 { 00245 mat temp(data.rows(), data.cols()); 00246 00247 for (int i = 0;i < temp.rows();i++) { 00248 for (int j = 0;j < temp.cols();j++) { 00249 temp(i, j) = std::arg(data(i, j)); 00250 } 00251 } 00252 00253 return temp; 00254 } 00255 00256 #ifdef _MSC_VER 00257 cvec conj(const cvec &x) 00258 { 00259 cvec temp(x.size()); 00260 00261 for (int i = 0; i < x.size(); i++) { 00262 temp(i) = std::conj(x(i)); 00263 } 00264 00265 return temp; 00266 } 00267 00268 cmat conj(const cmat &x) 00269 { 00270 cmat temp(x.rows(), x.cols()); 00271 00272 for (int i = 0; i < x.rows(); i++) { 00273 for (int j = 0; j < x.cols(); j++) { 00274 temp(i, j) = std::conj(x(i, j)); 00275 } 00276 } 00277 00278 return temp; 00279 } 00280 #endif 00281 00282 } // namespace itpp
Generated on Thu Apr 23 20:06:46 2009 for IT++ by Doxygen 1.5.8