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
Generated on Wed Jan 20 23:03:04 2010 for IT++ by Doxygen 1.6.2