00001 00033 #include <itpp/base/scalfunc.h> 00034 #include <itpp/base/vec.h> 00035 #include <cmath> 00036 #include <cstdlib> 00037 00038 00039 #ifndef HAVE_TGAMMA 00040 // "True" gamma function 00041 double tgamma(double x) 00042 { 00043 double s = (2.50662827510730 + 190.9551718930764 / (x + 1) 00044 - 216.8366818437280 / (x + 2) + 60.19441764023333 00045 / (x + 3) - 3.08751323928546 / (x + 4) + 0.00302963870525 00046 / (x + 5) - 0.00001352385959072596 / (x + 6)) / x; 00047 if (s < 0) 00048 return -exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(-s)); 00049 else 00050 return exp((x + 0.5) * log(x + 5.5) - x - 5.5 + log(s)); 00051 } 00052 #endif 00053 00054 #if !defined(HAVE_LGAMMA) || (HAVE_DECL_SIGNGAM != 1) 00055 // This global variable is normally declared in <cmath>, but not always 00056 int signgam; 00057 // Logarithm of an absolute value of gamma function 00058 double lgamma(double x) 00059 { 00060 double gam = tgamma(x); 00061 signgam = (gam < 0) ? -1 : 1; 00062 return log(fabs(gam)); 00063 } 00064 #endif 00065 00066 #ifndef HAVE_CBRT 00067 // Cubic root 00068 double cbrt(double x) 00069 { 00070 return pow(x, 1.0/3.0); 00071 } 00072 #endif 00073 00074 #ifndef HAVE_ASINH 00075 double asinh(double x) 00076 { 00077 return ((x>=0) ? log(x+sqrt(x*x+1)):-log(-x+sqrt(x*x+1))); 00078 } 00079 #endif 00080 00081 #ifndef HAVE_ACOSH 00082 double acosh(double x) 00083 { 00084 it_error_if(x<1,"acosh(x): x<1"); 00085 return log(x+sqrt(x*x-1)); 00086 } 00087 #endif 00088 00089 #ifndef HAVE_ATANH 00090 double atanh(double x) 00091 { 00092 it_error_if(fabs(x)>=1,"atanh(x): abs(x)>=1"); 00093 return 0.5*log((x+1)/(x-1)); 00094 } 00095 #endif 00096 00097 #ifdef _MSC_VER 00098 double erfc(double Y) 00099 { 00100 int ISW,I; 00101 double P[4],Q[3],P1[6],Q1[5],P2[4],Q2[3]; 00102 double XMIN,XLARGE,SQRPI; 00103 double X,RES,XSQ,XNUM,XDEN,XI,XBIG,ERFCret; 00104 P[1]=0.3166529; 00105 P[2]=1.722276; 00106 P[3]=21.38533; 00107 Q[1]=7.843746; 00108 Q[2]=18.95226; 00109 P1[1]=0.5631696; 00110 P1[2]=3.031799; 00111 P1[3]=6.865018; 00112 P1[4]=7.373888; 00113 P1[5]=4.318779e-5; 00114 Q1[1]=5.354217; 00115 Q1[2]=12.79553; 00116 Q1[3]=15.18491; 00117 Q1[4]=7.373961; 00118 P2[1]=5.168823e-2; 00119 P2[2]=0.1960690; 00120 P2[3]=4.257996e-2; 00121 Q2[1]=0.9214524; 00122 Q2[2]=0.1509421; 00123 XMIN=1.0E-5; 00124 XLARGE=4.1875E0; 00125 XBIG=9.0; 00126 SQRPI=0.5641896; 00127 X=Y; 00128 ISW=1; 00129 if (X<0) { 00130 ISW=-1; 00131 X=-X; 00132 } 00133 if (X<0.477) { 00134 if (X>=XMIN) { 00135 XSQ=X*X; 00136 XNUM=(P[1]*XSQ+P[2])*XSQ+P[3]; 00137 XDEN=(XSQ+Q[1])*XSQ+Q[2]; 00138 RES=X*XNUM/XDEN; 00139 } 00140 else RES=X*P[3]/Q[2]; 00141 if (ISW==-1) RES=-RES; 00142 RES=1.0-RES; 00143 goto slut; 00144 } 00145 if (X>4.0) { 00146 if (ISW>0) goto ulf; 00147 if (X<XLARGE) goto eva; 00148 RES=2.0; 00149 goto slut; 00150 } 00151 XSQ=X*X; 00152 XNUM=P1[5]*X+P1[1]; 00153 XDEN=X+Q1[1]; 00154 for(I=2;I<=4;I++) { 00155 XNUM=XNUM*X+P1[I]; 00156 XDEN=XDEN*X+Q1[I]; 00157 } 00158 RES=XNUM/XDEN; 00159 goto elin; 00160 ulf: if (X>XBIG) goto fred; 00161 eva: XSQ=X*X; 00162 XI=1.0/XSQ; 00163 XNUM=(P2[1]*XI+P2[2])*XI+P2[3]; 00164 XDEN=XI+Q2[1]*XI+Q2[2]; 00165 RES=(SQRPI+XI*XNUM/XDEN)/X; 00166 elin: RES=RES*exp(-XSQ); 00167 if (ISW==-1) RES=2.0-RES; 00168 goto slut; 00169 fred: RES=0.0; 00170 slut: ERFCret=RES; 00171 return ERFCret; 00172 } 00173 00174 double erf(double x) 00175 { 00176 return (1.0-erfc(x)); 00177 } 00178 #endif 00179 00180 namespace itpp { 00181 00182 double gamma(double x) 00183 { 00184 return tgamma(x); 00185 } 00186 00187 double Qfunc(double x) 00188 { 00189 return 0.5*erfc(x/1.41421356237310); 00190 } 00191 00192 double erfinv(double P) 00193 { 00194 double Y,A,B,X,Z,W,WI,SN,SD,F,Z2,SIGMA; 00195 double A1=-.5751703,A2=-1.896513,A3=-.5496261E-1; 00196 double B0=-.1137730,B1=-3.293474,B2=-2.374996,B3=-1.187515; 00197 double C0=-.1146666,C1=-.1314774,C2=-.2368201,C3=.5073975e-1; 00198 double D0=-44.27977,D1=21.98546,D2=-7.586103; 00199 double E0=-.5668422E-1,E1=.3937021,E2=-.3166501,E3=.6208963E-1; 00200 double F0=-6.266786,F1=4.666263,F2=-2.962883; 00201 double G0=.1851159E-3,G1=-.2028152E-2,G2=-.1498384,G3=.1078639E-1; 00202 double H0=.9952975E-1,H1=.5211733,H2=-.6888301E-1; 00203 // double RINFM=1.7014E+38; 00204 00205 X=P; 00206 SIGMA=sgn(X); 00207 it_error_if(X<-1 || X>1,"erfinv : argument out of bounds"); 00208 Z=fabs(X); 00209 if (Z>.85) { 00210 A=1-Z; 00211 B=Z; 00212 W=sqrt(-log(A+A*B)); 00213 if (W>=2.5) { 00214 if (W>=4.) { 00215 WI=1./W; 00216 SN=((G3*WI+G2)*WI+G1)*WI; 00217 SD=((WI+H2)*WI+H1)*WI+H0; 00218 F=W+W*(G0+SN/SD); 00219 } else { 00220 SN=((E3*W+E2)*W+E1)*W; 00221 SD=((W+F2)*W+F1)*W+F0; 00222 F=W+W*(E0+SN/SD); 00223 } 00224 } else { 00225 SN=((C3*W+C2)*W+C1)*W; 00226 SD=((W+D2)*W+D1)*W+D0; 00227 F=W+W*(C0+SN/SD); 00228 } 00229 } else { 00230 Z2=Z*Z; 00231 F=Z+Z*(B0+A1*Z2/(B1+Z2+A2/(B2+Z2+A3/(B3+Z2)))); 00232 } 00233 Y=SIGMA*F; 00234 return Y; 00235 } 00236 00237 00238 /* 00239 * Abramowitz and Stegun: Eq. (7.1.14) gives this continued fraction 00240 * for erfc(z) 00241 * 00242 * erfc(z) = sqrt(pi).exp(-z^2). 1 1/2 1 3/2 2 5/2 00243 * --- --- --- --- --- --- ... 00244 * z + z + z + z + z + z + 00245 * 00246 * This is evaluated using Lentz's method, as described in the 00247 * narative of Numerical Recipes in C. 00248 * 00249 * The continued fraction is true providing real(z) > 0. In practice 00250 * we like real(z) to be significantly greater than 0, say greater 00251 * than 0.5. 00252 */ 00253 std::complex<double> cerfc_continued_fraction(const std::complex<double>& z) 00254 { 00255 const double tiny = std::numeric_limits<double>::min(); 00256 00257 // first calculate z+ 1/2 1 00258 // --- --- ... 00259 // z + z + 00260 std::complex<double> f(z); 00261 std::complex<double> C(f); 00262 std::complex<double> D(0.0); 00263 std::complex<double> delta; 00264 double a; 00265 00266 a = 0.0; 00267 do { 00268 a += 0.5; 00269 D = z + a * D; 00270 C = z + a / C; 00271 if ((D.real() == 0.0) && (D.imag() == 0.0)) 00272 D = tiny; 00273 D = 1.0 / D; 00274 delta = C * D; 00275 f = f * delta; 00276 } while (abs(1.0 - delta) > eps); 00277 00278 // Do the first term of the continued fraction 00279 f = 1.0 / f; 00280 00281 // and do the final scaling 00282 f = f * exp(-z * z) / sqrt(pi); 00283 00284 return f; 00285 } 00286 00287 std::complex<double> cerf_continued_fraction(const std::complex<double>& z) 00288 { 00289 if (z.real() > 0) 00290 return 1.0 - cerfc_continued_fraction(z); 00291 else 00292 return -1.0 + cerfc_continued_fraction(-z); 00293 } 00294 00295 /* 00296 * Abramawitz and Stegun: Eq. (7.1.5) gives a series for erf(z) good 00297 * for all z, but converges faster for smallish abs(z), say abs(z) < 2. 00298 */ 00299 std::complex<double> cerf_series(const std::complex<double>& z) 00300 { 00301 const double tiny = std::numeric_limits<double>::min(); 00302 std::complex<double> sum(0.0); 00303 std::complex<double> term(z); 00304 std::complex<double> z2(z*z); 00305 00306 for (int n = 0; (n < 3) || (abs(term) > abs(sum) * tiny); n++) { 00307 sum += term / static_cast<double>(2 * n + 1); 00308 term *= -z2 / static_cast<double>(n + 1); 00309 } 00310 00311 return sum * 2.0 / sqrt(pi); 00312 } 00313 00314 /* 00315 * Numerical Recipes quotes a formula due to Rybicki for evaluating 00316 * Dawson's Integral: 00317 * 00318 * exp(-x^2) integral exp(t^2).dt = 1/sqrt(pi) lim sum exp(-(z-n.h)^2) / n 00319 * 0 to x h->0 n odd 00320 * 00321 * This can be adapted to erf(z). 00322 */ 00323 std::complex<double> cerf_rybicki(const std::complex<double>& z) 00324 { 00325 double h = 0.2; // numerical experiment suggests this is small enough 00326 00327 // choose an even n0, and then shift z->z-n0.h and n->n-h. 00328 // n0 is chosen so that real((z-n0.h)^2) is as small as possible. 00329 int n0 = 2 * static_cast<int>(z.imag() / (2 * h) + 0.5); 00330 00331 std::complex<double> z0(0.0, n0 * h); 00332 std::complex<double> zp(z - z0); 00333 std::complex<double> sum(0.0, 0.0); 00334 00335 // limits of sum chosen so that the end sums of the sum are 00336 // fairly small. In this case exp(-(35.h)^2)=5e-22 00337 for (int np = -35; np <= 35; np += 2) { 00338 std::complex<double> t(zp.real(), zp.imag() - np * h); 00339 std::complex<double> b(exp(t * t) / static_cast<double>(np + n0)); 00340 sum += b; 00341 } 00342 00343 sum *= 2.0 * exp(-z * z) / pi; 00344 00345 return std::complex<double>(-sum.imag(), sum.real()); 00346 } 00347 00348 /* 00349 * This function calculates a well known error function erf(z) for 00350 * complex z. Three methods are implemented. Which one is used 00351 * depends on z. 00352 */ 00353 std::complex<double> erf(const std::complex<double>& z) 00354 { 00355 // Use the method appropriate to size of z - 00356 // there probably ought to be an extra option for NaN z, or infinite z 00357 if (abs(z) < 2.0) 00358 return cerf_series(z); 00359 else { 00360 if (std::abs(z.real()) < 0.5) 00361 return cerf_rybicki(z); 00362 else 00363 return cerf_continued_fraction(z); 00364 } 00365 } 00366 00367 //Calculates factorial coefficient for index <= 170. 00368 double fact(int index) 00369 { 00370 it_error_if(index > 170,"\nThe function double factfp(int index) overflows if index > 170. \nUse your head instead!"); 00371 it_error_if(index < 0,"\nThe function double factfp(int index) cannot evaluate if index. < 0"); 00372 double prod = 1; 00373 for (int i=1; i<=index; i++) 00374 prod *= double(i); 00375 return prod; 00376 } 00377 00378 long mod(long k, long n) 00379 { 00380 if (n==0) { 00381 return k; 00382 } else { 00383 return (k - n * long(floor(double(k)/double(n))) ); 00384 } 00385 } 00386 00387 long gcd(long a, long b) 00388 { 00389 long v, u, t, q; 00390 00391 it_assert(a>=0,"long gcd(long a, long b): a and b must be non-negative integers"); 00392 it_assert(b>=0,"long gcd(long a, long b): a and b must be non-negative integers"); 00393 00394 u = std::abs(a); 00395 v = std::abs(b); 00396 while (v>0) { 00397 q = u / v; 00398 t = u - v*q; 00399 u = v; 00400 v = t; 00401 } 00402 return(u); 00403 } 00404 00405 00406 // Calculates binomial coefficient "n over k". 00407 double binom(int n, int k) 00408 { 00409 it_assert(k <= n, "binom(n, k): k can not be larger than n"); 00410 it_assert((n >= 0) && (k >= 0), "binom(n, k): n and k must be non-negative integers"); 00411 k = ((n - k) < k) ? n - k : k; 00412 00413 double out = n - k + 1; 00414 for (int i = 2; i <= k; i++) { 00415 out *= (i + n - k); 00416 out /= i; 00417 } 00418 return out; 00419 } 00420 00421 // Calculates binomial coefficient "n over k". 00422 int binom_i(int n, int k) 00423 { 00424 it_assert(k <= n, "binom_i(n, k): k can not be larger than n"); 00425 it_assert((n >= 0) && (k >= 0), "binom_i(n, k): n and k must be non-negative integers"); 00426 k = ((n - k) < k) ? n - k : k; 00427 00428 int out = n - k + 1; 00429 for (int i = 2; i <= k; i++) { 00430 out *= (i + n - k); 00431 out /= i; 00432 } 00433 return out; 00434 } 00435 00436 // Calculates the base 10-logarithm of the binomial coefficient "n over k". 00437 double log_binom(int n, int k) { 00438 it_assert(k <= n, "log_binom(n, k): k can not be larger than n"); 00439 it_assert((n >= 0) && (k >= 0), "log_binom(n, k): n and k must be non-negative integers"); 00440 k = ((n - k) < k) ? n - k : k; 00441 00442 double out = 0.0; 00443 for (int i = 1; i <= k; i++) 00444 out += log10(static_cast<double>(i + n - k)) 00445 - log10(static_cast<double>(i)); 00446 00447 return out; 00448 } 00449 00450 std::complex<double> round_to_zero(const std::complex<double>& x, 00451 double threshold) { 00452 return std::complex<double>(round_to_zero(x.real(), threshold), 00453 round_to_zero(x.imag(), threshold)); 00454 } 00455 00456 } // namespace itpp
Generated on Thu Aug 30 02:47:18 2007 for IT++ by Doxygen 1.5.3