IT++ Logo
error.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/math/error.h>
00030 #include <itpp/base/math/elem_math.h>
00031 #include <itpp/base/itcompat.h>
00032 
00033 
00034 namespace itpp
00035 {
00036 
00052 std::complex<double> cerfc_continued_fraction(const std::complex<double>& z)
00053 {
00054   const double tiny = std::numeric_limits<double>::min();
00055 
00056   // first calculate z+ 1/2   1
00057   //                    ---  --- ...
00058   //                    z +  z +
00059   std::complex<double> f(z);
00060   std::complex<double> C(f);
00061   std::complex<double> D(0.0);
00062   std::complex<double> delta;
00063   double a;
00064 
00065   a = 0.0;
00066   do {
00067     a += 0.5;
00068     D = z + a * D;
00069     C = z + a / C;
00070     if ((D.real() == 0.0) && (D.imag() == 0.0))
00071       D = tiny;
00072     D = 1.0 / D;
00073     delta = C * D;
00074     f = f * delta;
00075   }
00076   while (abs(1.0 - delta) > eps);
00077 
00078   // Do the first term of the continued fraction
00079   f = 1.0 / f;
00080 
00081   // and do the final scaling
00082   f = f * exp(-z * z) / std::sqrt(pi);
00083 
00084   return f;
00085 }
00086 
00088 std::complex<double> cerf_continued_fraction(const std::complex<double>& z)
00089 {
00090   if (z.real() > 0)
00091     return 1.0 - cerfc_continued_fraction(z);
00092   else
00093     return -1.0 + cerfc_continued_fraction(-z);
00094 }
00095 
00100 std::complex<double> cerf_series(const std::complex<double>& z)
00101 {
00102   const double tiny = std::numeric_limits<double>::min();
00103   std::complex<double> sum(0.0);
00104   std::complex<double> term(z);
00105   std::complex<double> z2(z*z);
00106 
00107   for (int n = 0; (n < 3) || (abs(term) > abs(sum) * tiny); n++) {
00108     sum += term / static_cast<double>(2 * n + 1);
00109     term *= -z2 / static_cast<double>(n + 1);
00110   }
00111 
00112   return sum * 2.0 / std::sqrt(pi);
00113 }
00114 
00124 std::complex<double> cerf_rybicki(const std::complex<double>& z)
00125 {
00126   double h = 0.2; // numerical experiment suggests this is small enough
00127 
00128   // choose an even n0, and then shift z->z-n0.h and n->n-h.
00129   // n0 is chosen so that real((z-n0.h)^2) is as small as possible.
00130   int n0 = 2 * static_cast<int>(z.imag() / (2 * h) + 0.5);
00131 
00132   std::complex<double> z0(0.0, n0 * h);
00133   std::complex<double> zp(z - z0);
00134   std::complex<double> sum(0.0, 0.0);
00135 
00136   // limits of sum chosen so that the end sums of the sum are
00137   // fairly small. In this case exp(-(35.h)^2)=5e-22
00138   for (int np = -35; np <= 35; np += 2) {
00139     std::complex<double> t(zp.real(), zp.imag() - np * h);
00140     std::complex<double> b(exp(t * t) / static_cast<double>(np + n0));
00141     sum += b;
00142   }
00143 
00144   sum *= 2.0 * exp(-z * z) / pi;
00145 
00146   return std::complex<double>(-sum.imag(), sum.real());
00147 }
00148 
00149 /*
00150  * This function calculates a well known error function erf(z) for
00151  * complex z. Three methods are implemented. Which one is used
00152  * depends on z.
00153  */
00154 std::complex<double> erf(const std::complex<double>& z)
00155 {
00156   // Use the method appropriate to size of z -
00157   // there probably ought to be an extra option for NaN z, or infinite z
00158   if (abs(z) < 2.0)
00159     return cerf_series(z);
00160   else {
00161     if (std::abs(z.real()) < 0.5)
00162       return cerf_rybicki(z);
00163     else
00164       return cerf_continued_fraction(z);
00165   }
00166 }
00167 
00168 
00169 double erfinv(double P)
00170 {
00171   double Y, A, B, X, Z, W, WI, SN, SD, F, Z2, SIGMA;
00172   double A1 = -.5751703, A2 = -1.896513, A3 = -.5496261E-1;
00173   double B0 = -.1137730, B1 = -3.293474, B2 = -2.374996, B3 = -1.187515;
00174   double C0 = -.1146666, C1 = -.1314774, C2 = -.2368201, C3 = .5073975e-1;
00175   double D0 = -44.27977, D1 = 21.98546, D2 = -7.586103;
00176   double E0 = -.5668422E-1, E1 = .3937021, E2 = -.3166501, E3 = .6208963E-1;
00177   double F0 = -6.266786, F1 = 4.666263, F2 = -2.962883;
00178   double G0 = .1851159E-3, G1 = -.2028152E-2, G2 = -.1498384, G3 = .1078639E-1;
00179   double H0 = .9952975E-1, H1 = .5211733, H2 = -.6888301E-1;
00180   // double RINFM=1.7014E+38;
00181 
00182   X = P;
00183   SIGMA = sign(X);
00184   it_error_if(X < -1 || X > 1, "erfinv : argument out of bounds");
00185   Z = fabs(X);
00186   if (Z > .85) {
00187     A = 1 - Z;
00188     B = Z;
00189     W = std::sqrt(-log(A + A * B));
00190     if (W >= 2.5) {
00191       if (W >= 4.) {
00192         WI = 1. / W;
00193         SN = ((G3 * WI + G2) * WI + G1) * WI;
00194         SD = ((WI + H2) * WI + H1) * WI + H0;
00195         F = W + W * (G0 + SN / SD);
00196       }
00197       else {
00198         SN = ((E3 * W + E2) * W + E1) * W;
00199         SD = ((W + F2) * W + F1) * W + F0;
00200         F = W + W * (E0 + SN / SD);
00201       }
00202     }
00203     else {
00204       SN = ((C3 * W + C2) * W + C1) * W;
00205       SD = ((W + D2) * W + D1) * W + D0;
00206       F = W + W * (C0 + SN / SD);
00207     }
00208   }
00209   else {
00210     Z2 = Z * Z;
00211     F = Z + Z * (B0 + A1 * Z2 / (B1 + Z2 + A2 / (B2 + Z2 + A3 / (B3 + Z2))));
00212   }
00213   Y = SIGMA * F;
00214   return Y;
00215 }
00216 
00217 double Qfunc(double x)
00218 {
00219   return (0.5 * ::erfc(x / 1.41421356237310));
00220 }
00221 
00222 
00223 // Error function
00224 vec erf(const vec &x) { return apply_function<double>(::erf, x); }
00225 mat erf(const mat &x) { return apply_function<double>(::erf, x); }
00226 cvec erf(const cvec &x)
00227 {
00228   return apply_function<std::complex<double> >(erf, x);
00229 }
00230 cmat erf(const cmat &x)
00231 {
00232   return apply_function<std::complex<double> >(erf, x);
00233 }
00234 
00235 // Inverse of error function
00236 vec erfinv(const vec &x) { return apply_function<double>(erfinv, x); }
00237 mat erfinv(const mat &x) { return apply_function<double>(erfinv, x); }
00238 
00239 // Complementary error function
00240 vec erfc(const vec &x) { return apply_function<double>(::erfc, x); }
00241 mat erfc(const mat &x) { return apply_function<double>(::erfc, x); }
00242 
00243 // Q-function
00244 vec Qfunc(const vec &x) { return apply_function<double>(Qfunc, x); }
00245 mat Qfunc(const mat &x) { return apply_function<double>(Qfunc, x); }
00246 
00247 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4