IT++ Logo Newcom Logo

scalfunc.cpp

Go to the documentation of this file.
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 = 1.0;
00414     for (int i = 1; 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 = 1;
00429     for (int i = 1; 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
SourceForge Logo

Generated on Sat Aug 25 23:40:25 2007 for IT++ by Doxygen 1.5.2