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:
00098   if (X > XBIG) goto fred;
00099 eva:
00100   XSQ = X * X;
00101   XI = 1.0 / XSQ;
00102   XNUM = (P2[1] * XI + P2[2]) * XI + P2[3];
00103   XDEN = XI + Q2[1] * XI + Q2[2];
00104   RES = (SQRPI + XI * XNUM / XDEN) / X;
00105 elin:
00106   RES = RES * exp(-XSQ);
00107   if (ISW == -1) RES = 2.0 - RES;
00108   goto slut;
00109 fred:
00110   RES = 0.0;
00111 slut:
00112   ERFCret = RES;
00113   return  ERFCret;
00114 }
00115 #endif
00116 
00117 #ifndef HAVE_ERF
00118 double erf(double x)
00119 {
00120   return (1.0 - ::erfc(x));
00121 }
00122 #endif
00123 
00124 
00125 namespace itpp
00126 {
00127 
00143 std::complex<double> cerfc_continued_fraction(const std::complex<double>& z)
00144 {
00145   const double tiny = std::numeric_limits<double>::min();
00146 
00147   // first calculate z+ 1/2   1
00148   //                    ---  --- ...
00149   //                    z +  z +
00150   std::complex<double> f(z);
00151   std::complex<double> C(f);
00152   std::complex<double> D(0.0);
00153   std::complex<double> delta;
00154   double a;
00155 
00156   a = 0.0;
00157   do {
00158     a += 0.5;
00159     D = z + a * D;
00160     C = z + a / C;
00161     if ((D.real() == 0.0) && (D.imag() == 0.0))
00162       D = tiny;
00163     D = 1.0 / D;
00164     delta = C * D;
00165     f = f * delta;
00166   }
00167   while (abs(1.0 - delta) > eps);
00168 
00169   // Do the first term of the continued fraction
00170   f = 1.0 / f;
00171 
00172   // and do the final scaling
00173   f = f * exp(-z * z) / std::sqrt(pi);
00174 
00175   return f;
00176 }
00177 
00179 std::complex<double> cerf_continued_fraction(const std::complex<double>& z)
00180 {
00181   if (z.real() > 0)
00182     return 1.0 - cerfc_continued_fraction(z);
00183   else
00184     return -1.0 + cerfc_continued_fraction(-z);
00185 }
00186 
00191 std::complex<double> cerf_series(const std::complex<double>& z)
00192 {
00193   const double tiny = std::numeric_limits<double>::min();
00194   std::complex<double> sum(0.0);
00195   std::complex<double> term(z);
00196   std::complex<double> z2(z*z);
00197 
00198   for (int n = 0; (n < 3) || (abs(term) > abs(sum) * tiny); n++) {
00199     sum += term / static_cast<double>(2 * n + 1);
00200     term *= -z2 / static_cast<double>(n + 1);
00201   }
00202 
00203   return sum * 2.0 / std::sqrt(pi);
00204 }
00205 
00215 std::complex<double> cerf_rybicki(const std::complex<double>& z)
00216 {
00217   double h = 0.2; // numerical experiment suggests this is small enough
00218 
00219   // choose an even n0, and then shift z->z-n0.h and n->n-h.
00220   // n0 is chosen so that real((z-n0.h)^2) is as small as possible.
00221   int n0 = 2 * static_cast<int>(z.imag() / (2 * h) + 0.5);
00222 
00223   std::complex<double> z0(0.0, n0 * h);
00224   std::complex<double> zp(z - z0);
00225   std::complex<double> sum(0.0, 0.0);
00226 
00227   // limits of sum chosen so that the end sums of the sum are
00228   // fairly small. In this case exp(-(35.h)^2)=5e-22
00229   for (int np = -35; np <= 35; np += 2) {
00230     std::complex<double> t(zp.real(), zp.imag() - np * h);
00231     std::complex<double> b(exp(t * t) / static_cast<double>(np + n0));
00232     sum += b;
00233   }
00234 
00235   sum *= 2.0 * exp(-z * z) / pi;
00236 
00237   return std::complex<double>(-sum.imag(), sum.real());
00238 }
00239 
00240 /*
00241  * This function calculates a well known error function erf(z) for
00242  * complex z. Three methods are implemented. Which one is used
00243  * depends on z.
00244  */
00245 std::complex<double> erf(const std::complex<double>& z)
00246 {
00247   // Use the method appropriate to size of z -
00248   // there probably ought to be an extra option for NaN z, or infinite z
00249   if (abs(z) < 2.0)
00250     return cerf_series(z);
00251   else {
00252     if (std::abs(z.real()) < 0.5)
00253       return cerf_rybicki(z);
00254     else
00255       return cerf_continued_fraction(z);
00256   }
00257 }
00258 
00259 
00260 double erfinv(double P)
00261 {
00262   double Y, A, B, X, Z, W, WI, SN, SD, F, Z2, SIGMA;
00263   double A1 = -.5751703, A2 = -1.896513, A3 = -.5496261E-1;
00264   double B0 = -.1137730, B1 = -3.293474, B2 = -2.374996, B3 = -1.187515;
00265   double C0 = -.1146666, C1 = -.1314774, C2 = -.2368201, C3 = .5073975e-1;
00266   double D0 = -44.27977, D1 = 21.98546, D2 = -7.586103;
00267   double E0 = -.5668422E-1, E1 = .3937021, E2 = -.3166501, E3 = .6208963E-1;
00268   double F0 = -6.266786, F1 = 4.666263, F2 = -2.962883;
00269   double G0 = .1851159E-3, G1 = -.2028152E-2, G2 = -.1498384, G3 = .1078639E-1;
00270   double H0 = .9952975E-1, H1 = .5211733, H2 = -.6888301E-1;
00271   // double RINFM=1.7014E+38;
00272 
00273   X = P;
00274   SIGMA = sign(X);
00275   it_error_if(X < -1 || X > 1, "erfinv : argument out of bounds");
00276   Z = fabs(X);
00277   if (Z > .85) {
00278     A = 1 - Z;
00279     B = Z;
00280     W = std::sqrt(-log(A + A * B));
00281     if (W >= 2.5) {
00282       if (W >= 4.) {
00283         WI = 1. / W;
00284         SN = ((G3 * WI + G2) * WI + G1) * WI;
00285         SD = ((WI + H2) * WI + H1) * WI + H0;
00286         F = W + W * (G0 + SN / SD);
00287       }
00288       else {
00289         SN = ((E3 * W + E2) * W + E1) * W;
00290         SD = ((W + F2) * W + F1) * W + F0;
00291         F = W + W * (E0 + SN / SD);
00292       }
00293     }
00294     else {
00295       SN = ((C3 * W + C2) * W + C1) * W;
00296       SD = ((W + D2) * W + D1) * W + D0;
00297       F = W + W * (C0 + SN / SD);
00298     }
00299   }
00300   else {
00301     Z2 = Z * Z;
00302     F = Z + Z * (B0 + A1 * Z2 / (B1 + Z2 + A2 / (B2 + Z2 + A3 / (B3 + Z2))));
00303   }
00304   Y = SIGMA * F;
00305   return Y;
00306 }
00307 
00308 double Qfunc(double x)
00309 {
00310   return (0.5 * ::erfc(x / 1.41421356237310));
00311 }
00312 
00313 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Tue Feb 2 09:33:29 2010 for IT++ by Doxygen 1.6.2