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
Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4