IT++ Logo

error.cpp

Go to the documentation of this file.
00001 
00030 #include <itpp/base/math/error.h>
00031 #include <itpp/base/math/elem_math.h>
00032 
00033 
00034 #ifndef HAVE_ERFC
00035 double erfc(double Y)
00036 {
00037   int  ISW,I;
00038   double P[4],Q[3],P1[6],Q1[5],P2[4],Q2[3];
00039   double XMIN,XLARGE,SQRPI;
00040   double X,RES,XSQ,XNUM,XDEN,XI,XBIG,ERFCret;
00041   P[1]=0.3166529;
00042   P[2]=1.722276;
00043   P[3]=21.38533;
00044   Q[1]=7.843746;
00045   Q[2]=18.95226;
00046   P1[1]=0.5631696;
00047   P1[2]=3.031799;
00048   P1[3]=6.865018;
00049   P1[4]=7.373888;
00050   P1[5]=4.318779e-5;
00051   Q1[1]=5.354217;
00052   Q1[2]=12.79553;
00053   Q1[3]=15.18491;
00054   Q1[4]=7.373961;
00055   P2[1]=5.168823e-2;
00056   P2[2]=0.1960690;
00057   P2[3]=4.257996e-2;
00058   Q2[1]=0.9214524;
00059   Q2[2]=0.1509421;
00060   XMIN=1.0E-5;
00061   XLARGE=4.1875E0;
00062   XBIG=9.0;
00063   SQRPI=0.5641896;
00064   X=Y;
00065   ISW=1;
00066   if (X<0) {
00067     ISW=-1;
00068     X=-X;
00069   }
00070   if (X<0.477) {
00071     if (X>=XMIN) {
00072       XSQ=X*X;
00073       XNUM=(P[1]*XSQ+P[2])*XSQ+P[3];
00074       XDEN=(XSQ+Q[1])*XSQ+Q[2];
00075       RES=X*XNUM/XDEN;
00076     }
00077     else RES=X*P[3]/Q[2];
00078     if (ISW==-1) RES=-RES;
00079     RES=1.0-RES;
00080     goto slut;
00081   }
00082   if (X>4.0) {
00083     if (ISW>0) goto ulf;
00084     if (X<XLARGE) goto eva;
00085     RES=2.0;
00086     goto slut;
00087   }
00088   XSQ=X*X;
00089   XNUM=P1[5]*X+P1[1];
00090   XDEN=X+Q1[1];
00091   for(I=2;I<=4;I++) {
00092     XNUM=XNUM*X+P1[I];
00093     XDEN=XDEN*X+Q1[I];
00094   }
00095   RES=XNUM/XDEN;
00096   goto elin;
00097  ulf:   if (X>XBIG) goto fred;
00098  eva:   XSQ=X*X;
00099   XI=1.0/XSQ;
00100   XNUM=(P2[1]*XI+P2[2])*XI+P2[3];
00101   XDEN=XI+Q2[1]*XI+Q2[2];
00102   RES=(SQRPI+XI*XNUM/XDEN)/X;
00103  elin:  RES=RES*exp(-XSQ);
00104   if (ISW==-1) RES=2.0-RES;
00105   goto slut;
00106  fred:  RES=0.0;
00107  slut:  ERFCret=RES;
00108   return  ERFCret;
00109 }
00110 #endif
00111 
00112 #ifndef HAVE_ERF
00113 double erf(double x)
00114 {
00115   return (1.0 - ::erfc(x));
00116 }
00117 #endif
00118 
00119 
00120 namespace itpp {
00121 
00137   std::complex<double> cerfc_continued_fraction(const std::complex<double>& z)
00138   {
00139     const double tiny = std::numeric_limits<double>::min();
00140 
00141     // first calculate z+ 1/2   1
00142     //                    ---  --- ...
00143     //                    z +  z +
00144     std::complex<double> f(z);
00145     std::complex<double> C(f);
00146     std::complex<double> D(0.0);
00147     std::complex<double> delta;
00148     double a;
00149 
00150     a = 0.0;
00151     do {
00152       a += 0.5;
00153       D = z + a * D;
00154       C = z + a / C;
00155       if ((D.real() == 0.0) && (D.imag() == 0.0))
00156         D = tiny;
00157       D = 1.0 / D;
00158       delta = C * D;
00159       f = f * delta;
00160     } while (abs(1.0 - delta) > eps);
00161 
00162     // Do the first term of the continued fraction
00163     f = 1.0 / f;
00164 
00165     // and do the final scaling
00166     f = f * exp(-z * z) / std::sqrt(pi);
00167 
00168     return f;
00169   }
00170 
00172   std::complex<double> cerf_continued_fraction(const std::complex<double>& z)
00173   {
00174     if (z.real() > 0)
00175       return 1.0 - cerfc_continued_fraction(z);
00176     else
00177       return -1.0 + cerfc_continued_fraction(-z);
00178   }
00179 
00184   std::complex<double> cerf_series(const std::complex<double>& z)
00185   {
00186     const double tiny = std::numeric_limits<double>::min();
00187     std::complex<double> sum(0.0);
00188     std::complex<double> term(z);
00189     std::complex<double> z2(z*z);
00190 
00191     for (int n = 0; (n < 3) || (abs(term) > abs(sum) * tiny); n++) {
00192       sum += term / static_cast<double>(2 * n + 1);
00193       term *= -z2 / static_cast<double>(n + 1);
00194     }
00195 
00196     return sum * 2.0 / std::sqrt(pi);
00197   }
00198 
00208   std::complex<double> cerf_rybicki(const std::complex<double>& z)
00209   {
00210     double h = 0.2; // numerical experiment suggests this is small enough
00211 
00212     // choose an even n0, and then shift z->z-n0.h and n->n-h.
00213     // n0 is chosen so that real((z-n0.h)^2) is as small as possible.
00214     int n0 = 2 * static_cast<int>(z.imag() / (2 * h) + 0.5);
00215 
00216     std::complex<double> z0(0.0, n0 * h);
00217     std::complex<double> zp(z - z0);
00218     std::complex<double> sum(0.0, 0.0);
00219 
00220     // limits of sum chosen so that the end sums of the sum are
00221     // fairly small. In this case exp(-(35.h)^2)=5e-22
00222     for (int np = -35; np <= 35; np += 2) {
00223       std::complex<double> t(zp.real(), zp.imag() - np * h);
00224       std::complex<double> b(exp(t * t) / static_cast<double>(np + n0));
00225       sum += b;
00226     }
00227 
00228     sum *= 2.0 * exp(-z * z) / pi;
00229 
00230     return std::complex<double>(-sum.imag(), sum.real());
00231   }
00232 
00233   /*
00234    * This function calculates a well known error function erf(z) for
00235    * complex z. Three methods are implemented. Which one is used
00236    * depends on z.
00237    */
00238   std::complex<double> erf(const std::complex<double>& z)
00239   {
00240     // Use the method appropriate to size of z -
00241     // there probably ought to be an extra option for NaN z, or infinite z
00242     if (abs(z) < 2.0)
00243       return cerf_series(z);
00244     else {
00245       if (std::abs(z.real()) < 0.5)
00246         return cerf_rybicki(z);
00247       else
00248         return cerf_continued_fraction(z);
00249     }
00250   }
00251 
00252 
00253   double erfinv(double P)
00254   {
00255     double  Y,A,B,X,Z,W,WI,SN,SD,F,Z2,SIGMA;
00256     double  A1=-.5751703,A2=-1.896513,A3=-.5496261E-1;
00257     double  B0=-.1137730,B1=-3.293474,B2=-2.374996,B3=-1.187515;
00258     double  C0=-.1146666,C1=-.1314774,C2=-.2368201,C3=.5073975e-1;
00259     double  D0=-44.27977,D1=21.98546,D2=-7.586103;
00260     double  E0=-.5668422E-1,E1=.3937021,E2=-.3166501,E3=.6208963E-1;
00261     double  F0=-6.266786,F1=4.666263,F2=-2.962883;
00262     double  G0=.1851159E-3,G1=-.2028152E-2,G2=-.1498384,G3=.1078639E-1;
00263     double  H0=.9952975E-1,H1=.5211733,H2=-.6888301E-1;
00264     //  double  RINFM=1.7014E+38;
00265 
00266     X=P;
00267     SIGMA = sign(X);
00268     it_error_if(X<-1 || X>1,"erfinv : argument out of bounds");
00269     Z=fabs(X);
00270     if (Z>.85) {
00271       A=1-Z;
00272       B=Z;
00273       W = std::sqrt(-log(A+A*B));
00274       if (W>=2.5) {
00275   if (W>=4.) {
00276     WI=1./W;
00277     SN=((G3*WI+G2)*WI+G1)*WI;
00278     SD=((WI+H2)*WI+H1)*WI+H0;
00279     F=W+W*(G0+SN/SD);
00280   } else {
00281     SN=((E3*W+E2)*W+E1)*W;
00282     SD=((W+F2)*W+F1)*W+F0;
00283     F=W+W*(E0+SN/SD);
00284   }
00285       } else {
00286   SN=((C3*W+C2)*W+C1)*W;
00287   SD=((W+D2)*W+D1)*W+D0;
00288   F=W+W*(C0+SN/SD);
00289       }
00290     } else {
00291       Z2=Z*Z;
00292       F=Z+Z*(B0+A1*Z2/(B1+Z2+A2/(B2+Z2+A3/(B3+Z2))));
00293     }
00294     Y=SIGMA*F;
00295     return Y;
00296   }
00297 
00298   double Qfunc(double x)
00299   {
00300     return (0.5 * ::erfc(x / 1.41421356237310));
00301   }
00302 
00303 } // namespace itpp
SourceForge Logo

Generated on Sun Sep 14 18:54:52 2008 for IT++ by Doxygen 1.5.6