IT++ Logo
itcompat.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/itcompat.h>
00030 #include <itpp/base/itassert.h>
00031 #include <limits>
00032 
00034 
00035 #ifndef HAVE_TGAMMA
00036 // "True" gamma function
00037 double tgamma(double x)
00038 {
00039   double s = (2.50662827510730 + 190.9551718930764 / (x + 1)
00040               - 216.8366818437280 / (x + 2) + 60.19441764023333
00041               / (x + 3) - 3.08751323928546 / (x + 4) + 0.00302963870525
00042               / (x + 5) - 0.00001352385959072596 / (x + 6)) / x;
00043   if (s < 0)
00044     return -exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(-s));
00045   else
00046     return exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(s));
00047 }
00048 #endif
00049 
00050 #if !defined(HAVE_LGAMMA) || (HAVE_DECL_SIGNGAM != 1)
00051 // The sign of the Gamma function is returned in the external integer
00052 // signgam declared in <math.h>. It is 1 when the Gamma function is positive
00053 // or zero, -1 when it is negative. However, MinGW definition of lgamma()
00054 // function does not use the global signgam variable.
00055 int signgam;
00056 // Logarithm of an absolute value of gamma function
00057 double lgamma(double x)
00058 {
00059   double gam = tgamma(x);
00060   signgam = (gam < 0) ? -1 : 1;
00061   return log(fabs(gam));
00062 }
00063 #endif
00064 
00065 #ifndef HAVE_CBRT
00066 // Cubic root
00067 double cbrt(double x) { return std::pow(x, 1.0 / 3.0); }
00068 #endif
00069 
00070 
00071 #ifndef HAVE_EXPM1
00072 #ifndef M_LN2
00073 #define M_LN2 0.69314718055994530941723212146 // ln(2)
00074 #endif
00075 // Implementation taken from GSL (http://www.gnu.org/software/gsl/)
00076 double expm1(double x)
00077 {
00078   /* FIXME: this should be improved */
00079   if (std::fabs(x) < M_LN2) {
00080     /* Compute the taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ... */
00081     double i = 1.0;
00082     double sum = x;
00083     double term = x / 1.0;
00084     do {
00085       i++;
00086       term *= x / i;
00087       sum += term;
00088     }
00089     while (std::fabs(term) > (std::fabs(sum)
00090                               * std::numeric_limits<double>::epsilon()));
00091     return sum;
00092   }
00093   else {
00094     return std::exp(x) - 1.0;
00095   }
00096 }
00097 #endif // HAVE_EXPM1
00098 
00099 
00100 #ifndef HAVE_ERFC
00101 // Complementary error function
00102 double erfc(double Y)
00103 {
00104   int  ISW, I;
00105   double P[4], Q[3], P1[6], Q1[5], P2[4], Q2[3];
00106   double XMIN, XLARGE, SQRPI;
00107   double X, RES, XSQ, XNUM, XDEN, XI, XBIG, ERFCret;
00108   P[1] = 0.3166529;
00109   P[2] = 1.722276;
00110   P[3] = 21.38533;
00111   Q[1] = 7.843746;
00112   Q[2] = 18.95226;
00113   P1[1] = 0.5631696;
00114   P1[2] = 3.031799;
00115   P1[3] = 6.865018;
00116   P1[4] = 7.373888;
00117   P1[5] = 4.318779e-5;
00118   Q1[1] = 5.354217;
00119   Q1[2] = 12.79553;
00120   Q1[3] = 15.18491;
00121   Q1[4] = 7.373961;
00122   P2[1] = 5.168823e-2;
00123   P2[2] = 0.1960690;
00124   P2[3] = 4.257996e-2;
00125   Q2[1] = 0.9214524;
00126   Q2[2] = 0.1509421;
00127   XMIN = 1.0E-5;
00128   XLARGE = 4.1875E0;
00129   XBIG = 9.0;
00130   SQRPI = 0.5641896;
00131   X = Y;
00132   ISW = 1;
00133   if (X < 0) {
00134     ISW = -1;
00135     X = -X;
00136   }
00137   if (X < 0.477) {
00138     if (X >= XMIN) {
00139       XSQ = X * X;
00140       XNUM = (P[1] * XSQ + P[2]) * XSQ + P[3];
00141       XDEN = (XSQ + Q[1]) * XSQ + Q[2];
00142       RES = X * XNUM / XDEN;
00143     }
00144     else RES = X * P[3] / Q[2];
00145     if (ISW == -1) RES = -RES;
00146     RES = 1.0 - RES;
00147     goto slut;
00148   }
00149   if (X > 4.0) {
00150     if (ISW > 0) goto ulf;
00151     if (X < XLARGE) goto eva;
00152     RES = 2.0;
00153     goto slut;
00154   }
00155   XSQ = X * X;
00156   XNUM = P1[5] * X + P1[1];
00157   XDEN = X + Q1[1];
00158   for (I = 2;I <= 4;I++) {
00159     XNUM = XNUM * X + P1[I];
00160     XDEN = XDEN * X + Q1[I];
00161   }
00162   RES = XNUM / XDEN;
00163   goto elin;
00164 ulf:
00165   if (X > XBIG) goto fred;
00166 eva:
00167   XSQ = X * X;
00168   XI = 1.0 / XSQ;
00169   XNUM = (P2[1] * XI + P2[2]) * XI + P2[3];
00170   XDEN = XI + Q2[1] * XI + Q2[2];
00171   RES = (SQRPI + XI * XNUM / XDEN) / X;
00172 elin:
00173   RES = RES * exp(-XSQ);
00174   if (ISW == -1) RES = 2.0 - RES;
00175   goto slut;
00176 fred:
00177   RES = 0.0;
00178 slut:
00179   ERFCret = RES;
00180   return ERFCret;
00181 }
00182 #endif
00183 
00184 
00185 #ifndef HAVE_ASINH
00186 // Arcus sinhyp
00187 double asinh(double x)
00188 {
00189   return ((x >= 0) ? std::log(x + std::sqrt(x * x + 1))
00190           : -std::log(-x + std::sqrt(x * x + 1)));
00191 }
00192 #endif
00193 
00194 #ifndef HAVE_ACOSH
00195 // Arcus coshyp
00196 double acosh(double x)
00197 {
00198   it_error_if(x < 1, "acosh(): Argument out of range");
00199   return std::log(x + std::sqrt(x * x - 1.0));
00200 }
00201 #endif
00202 
00203 #ifndef HAVE_ATANH
00204 // Arcus tanhyp
00205 double atanh(double x)
00206 {
00207   it_error_if(std::fabs(x) >= 1, "atanh(): Argument out of range");
00208   return 0.5 * std::log((x + 1) / (x - 1));
00209 }
00210 #endif
00211 
00212 
00213 #ifndef HAVE_RINT
00214 double rint(double x)
00215 {
00216   // zero or NaN case
00217   if ((x == 0.0) || (x != x))
00218     return x;
00219 
00220   // negative case
00221   bool neg = false;
00222   if (x < 0.0) {
00223     x = -x;
00224     neg = true;
00225   }
00226 
00227   double y = std::floor(x + 0.5);
00228   int i = static_cast<int>(y);
00229   if ((y - x >= 0.5) && (i & 1))
00230     --y;
00231 
00232   return neg ? -y : y;
00233 }
00234 #endif // HAVE_RINT
00235 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4