IT++ Logo

elem_math.cpp

Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Feb 9 2011 13:47:16 for IT++ by Doxygen 1.7.3