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   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
SourceForge Logo

Generated on Sun Sep 14 18:52:34 2008 for IT++ by Doxygen 1.5.6