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
Generated on Sun Sep 14 18:54:52 2008 for IT++ by Doxygen 1.5.6