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
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4