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