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 vec sqr(const cvec &data) 00072 { 00073 vec temp(data.length()); 00074 for (int i = 0; i < data.length(); i++) 00075 temp(i) = sqr(data(i)); 00076 return temp; 00077 } 00078 00079 mat sqr(const cmat &data) 00080 { 00081 mat temp(data.rows(), data.cols()); 00082 for (int i = 0; i < temp.rows(); i++) { 00083 for (int j = 0; j < temp.cols(); j++) { 00084 temp(i, j) = sqr(data(i, j)); 00085 } 00086 } 00087 return temp; 00088 } 00089 00090 vec abs(const cvec &data) 00091 { 00092 vec temp(data.length()); 00093 00094 for (int i=0;i<data.length();i++) 00095 temp[i]=std::abs(data[i]); 00096 00097 return temp; 00098 } 00099 00100 mat abs(const cmat &data) 00101 { 00102 mat temp(data.rows(),data.cols()); 00103 00104 for (int i=0;i<temp.rows();i++) { 00105 for (int j=0;j<temp.cols();j++) { 00106 temp(i,j)=std::abs(data(i,j)); 00107 } 00108 } 00109 00110 return temp; 00111 } 00112 00113 // Calculates factorial coefficient for index <= 170. 00114 double fact(int index) 00115 { 00116 it_error_if(index > 170, "fact(int index): Function overflows if index > 170."); 00117 it_error_if(index < 0, "fact(int index): index must be non-negative integer"); 00118 double prod = 1; 00119 for (int i = 1; i <= index; i++) 00120 prod *= static_cast<double>(i); 00121 return prod; 00122 } 00123 00124 // Calculates binomial coefficient "n over k". 00125 double binom(int n, int k) 00126 { 00127 it_assert(k <= n, "binom(n, k): k can not be larger than n"); 00128 it_assert((n >= 0) && (k >= 0), "binom(n, k): n and k must be non-negative integers"); 00129 k = ((n - k) < k) ? n - k : k; 00130 00131 double out = 1.0; 00132 for (int i = 1; i <= k; ++i) { 00133 out *= (i + n - k); 00134 out /= i; 00135 } 00136 return out; 00137 } 00138 00139 // Calculates binomial coefficient "n over k". 00140 int binom_i(int n, int k) 00141 { 00142 it_assert(k <= n, "binom_i(n, k): k can not be larger than n"); 00143 it_assert((n >= 0) && (k >= 0), "binom_i(n, k): n and k must be non-negative integers"); 00144 k = ((n - k) < k) ? n - k : k; 00145 00146 int out = 1; 00147 for (int i = 1; i <= k; ++i) { 00148 out *= (i + n - k); 00149 out /= i; 00150 } 00151 return out; 00152 } 00153 00154 // Calculates the base 10-logarithm of the binomial coefficient "n over k". 00155 double log_binom(int n, int k) { 00156 it_assert(k <= n, "log_binom(n, k): k can not be larger than n"); 00157 it_assert((n >= 0) && (k >= 0), "log_binom(n, k): n and k must be non-negative integers"); 00158 k = ((n - k) < k) ? n - k : k; 00159 00160 double out = 0.0; 00161 for (int i = 1; i <= k; i++) 00162 out += log10(static_cast<double>(i + n - k)) 00163 - log10(static_cast<double>(i)); 00164 00165 return out; 00166 } 00167 00168 // Calculates the greatest common divisor 00169 int gcd(int a, int b) 00170 { 00171 it_assert((a >= 0) && (b >= 0),"gcd(a, b): a and b must be non-negative integers"); 00172 int v, u, t, q; 00173 00174 u = a; 00175 v = b; 00176 while (v > 0) { 00177 q = u / v; 00178 t = u - v*q; 00179 u = v; 00180 v = t; 00181 } 00182 return u; 00183 } 00184 00185 00186 vec real(const cvec &data) 00187 { 00188 vec temp(data.length()); 00189 00190 for (int i=0;i<data.length();i++) 00191 temp[i]=data[i].real(); 00192 00193 return temp; 00194 } 00195 00196 mat real(const cmat &data) 00197 { 00198 mat temp(data.rows(),data.cols()); 00199 00200 for (int i=0;i<temp.rows();i++) { 00201 for (int j=0;j<temp.cols();j++) { 00202 temp(i,j)=data(i,j).real(); 00203 } 00204 } 00205 00206 return temp; 00207 } 00208 00209 vec imag(const cvec &data) 00210 { 00211 vec temp(data.length()); 00212 00213 for (int i=0;i<data.length();i++) 00214 temp[i]=data[i].imag(); 00215 return temp; 00216 } 00217 00218 mat imag(const cmat &data) 00219 { 00220 mat temp(data.rows(),data.cols()); 00221 00222 for (int i=0;i<temp.rows();i++) { 00223 for (int j=0;j<temp.cols();j++) { 00224 temp(i,j)=data(i,j).imag(); 00225 } 00226 } 00227 00228 return temp; 00229 } 00230 00231 vec arg(const cvec &data) 00232 { 00233 vec temp(data.length()); 00234 00235 for (int i=0;i<data.length();i++) 00236 temp[i]=std::arg(data[i]); 00237 00238 return temp; 00239 } 00240 00241 mat arg(const cmat &data) 00242 { 00243 mat temp(data.rows(),data.cols()); 00244 00245 for (int i=0;i<temp.rows();i++) { 00246 for (int j=0;j<temp.cols();j++) { 00247 temp(i,j)=std::arg(data(i,j)); 00248 } 00249 } 00250 00251 return temp; 00252 } 00253 00254 #ifdef _MSC_VER 00255 cvec conj(const cvec &x) { 00256 cvec temp(x.size()); 00257 00258 for (int i = 0; i < x.size(); i++) { 00259 temp(i) = std::conj(x(i)); 00260 } 00261 00262 return temp; 00263 } 00264 00265 cmat conj(const cmat &x) { 00266 cmat temp(x.rows(), x.cols()); 00267 00268 for (int i = 0; i < x.rows(); i++) { 00269 for (int j = 0; j < x.cols(); j++) { 00270 temp(i, j) = std::conj(x(i, j)); 00271 } 00272 } 00273 00274 return temp; 00275 } 00276 #endif 00277 00278 } // namespace itpp
Generated on Sun Dec 9 17:26:17 2007 for IT++ by Doxygen 1.5.4