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 
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
SourceForge Logo

Generated on Thu Apr 19 14:14:57 2007 for IT++ by Doxygen 1.5.1