IT++ Logo
error.cpp
Go to the documentation of this file.
1 
29 #include <itpp/base/math/error.h>
31 #include <itpp/base/itcompat.h>
32 
33 
34 namespace itpp
35 {
36 
52 std::complex<double> cerfc_continued_fraction(const std::complex<double>& z)
53 {
54  const double tiny = std::numeric_limits<double>::min();
55 
56  // first calculate z+ 1/2 1
57  // --- --- ...
58  // z + z +
59  std::complex<double> f(z);
60  std::complex<double> C(f);
61  std::complex<double> D(0.0);
62  std::complex<double> delta;
63  double a;
64 
65  a = 0.0;
66  do {
67  a += 0.5;
68  D = z + a * D;
69  C = z + a / C;
70  if ((D.real() == 0.0) && (D.imag() == 0.0))
71  D = tiny;
72  D = 1.0 / D;
73  delta = C * D;
74  f = f * delta;
75  }
76  while (abs(1.0 - delta) > eps);
77 
78  // Do the first term of the continued fraction
79  f = 1.0 / f;
80 
81  // and do the final scaling
82  f = f * exp(-z * z) / std::sqrt(pi);
83 
84  return f;
85 }
86 
88 std::complex<double> cerf_continued_fraction(const std::complex<double>& z)
89 {
90  if (z.real() > 0)
91  return 1.0 - cerfc_continued_fraction(z);
92  else
93  return -1.0 + cerfc_continued_fraction(-z);
94 }
95 
100 std::complex<double> cerf_series(const std::complex<double>& z)
101 {
102  const double tiny = std::numeric_limits<double>::min();
103  std::complex<double> sum(0.0);
104  std::complex<double> term(z);
105  std::complex<double> z2(z*z);
106 
107  for (int n = 0; (n < 3) || (abs(term) > abs(sum) * tiny); n++) {
108  sum += term / static_cast<double>(2 * n + 1);
109  term *= -z2 / static_cast<double>(n + 1);
110  }
111 
112  return sum * 2.0 / std::sqrt(pi);
113 }
114 
124 std::complex<double> cerf_rybicki(const std::complex<double>& z)
125 {
126  double h = 0.2; // numerical experiment suggests this is small enough
127 
128  // choose an even n0, and then shift z->z-n0.h and n->n-h.
129  // n0 is chosen so that real((z-n0.h)^2) is as small as possible.
130  int n0 = 2 * static_cast<int>(z.imag() / (2 * h) + 0.5);
131 
132  std::complex<double> z0(0.0, n0 * h);
133  std::complex<double> zp(z - z0);
134  std::complex<double> sum(0.0, 0.0);
135 
136  // limits of sum chosen so that the end sums of the sum are
137  // fairly small. In this case exp(-(35.h)^2)=5e-22
138  for (int np = -35; np <= 35; np += 2) {
139  std::complex<double> t(zp.real(), zp.imag() - np * h);
140  std::complex<double> b(exp(t * t) / static_cast<double>(np + n0));
141  sum += b;
142  }
143 
144  sum *= 2.0 * exp(-z * z) / pi;
145 
146  return std::complex<double>(-sum.imag(), sum.real());
147 }
148 
149 /*
150  * This function calculates a well known error function erf(z) for
151  * complex z. Three methods are implemented. Which one is used
152  * depends on z.
153  */
154 std::complex<double> erf(const std::complex<double>& z)
155 {
156  // Use the method appropriate to size of z -
157  // there probably ought to be an extra option for NaN z, or infinite z
158  if (abs(z) < 2.0)
159  return cerf_series(z);
160  else {
161  if (std::abs(z.real()) < 0.5)
162  return cerf_rybicki(z);
163  else
164  return cerf_continued_fraction(z);
165  }
166 }
167 
168 
169 double erfinv(double P)
170 {
171  double Y, A, B, X, Z, W, WI, SN, SD, F, Z2, SIGMA;
172  double A1 = -.5751703, A2 = -1.896513, A3 = -.5496261E-1;
173  double B0 = -.1137730, B1 = -3.293474, B2 = -2.374996, B3 = -1.187515;
174  double C0 = -.1146666, C1 = -.1314774, C2 = -.2368201, C3 = .5073975e-1;
175  double D0 = -44.27977, D1 = 21.98546, D2 = -7.586103;
176  double E0 = -.5668422E-1, E1 = .3937021, E2 = -.3166501, E3 = .6208963E-1;
177  double F0 = -6.266786, F1 = 4.666263, F2 = -2.962883;
178  double G0 = .1851159E-3, G1 = -.2028152E-2, G2 = -.1498384, G3 = .1078639E-1;
179  double H0 = .9952975E-1, H1 = .5211733, H2 = -.6888301E-1;
180  // double RINFM=1.7014E+38;
181 
182  X = P;
183  SIGMA = sign(X);
184  it_error_if(X < -1 || X > 1, "erfinv : argument out of bounds");
185  Z = fabs(X);
186  if (Z > .85) {
187  A = 1 - Z;
188  B = Z;
189  W = std::sqrt(-log(A + A * B));
190  if (W >= 2.5) {
191  if (W >= 4.) {
192  WI = 1. / W;
193  SN = ((G3 * WI + G2) * WI + G1) * WI;
194  SD = ((WI + H2) * WI + H1) * WI + H0;
195  F = W + W * (G0 + SN / SD);
196  }
197  else {
198  SN = ((E3 * W + E2) * W + E1) * W;
199  SD = ((W + F2) * W + F1) * W + F0;
200  F = W + W * (E0 + SN / SD);
201  }
202  }
203  else {
204  SN = ((C3 * W + C2) * W + C1) * W;
205  SD = ((W + D2) * W + D1) * W + D0;
206  F = W + W * (C0 + SN / SD);
207  }
208  }
209  else {
210  Z2 = Z * Z;
211  F = Z + Z * (B0 + A1 * Z2 / (B1 + Z2 + A2 / (B2 + Z2 + A3 / (B3 + Z2))));
212  }
213  Y = SIGMA * F;
214  return Y;
215 }
216 
217 double Qfunc(double x)
218 {
219  return (0.5 * ::erfc(x / 1.41421356237310));
220 }
221 
222 
223 // Error function
224 vec erf(const vec &x) { return apply_function<double>(::erf, x); }
225 mat erf(const mat &x) { return apply_function<double>(::erf, x); }
226 cvec erf(const cvec &x)
227 {
228  return apply_function<std::complex<double> >(erf, x);
229 }
230 cmat erf(const cmat &x)
231 {
232  return apply_function<std::complex<double> >(erf, x);
233 }
234 
235 // Inverse of error function
236 vec erfinv(const vec &x) { return apply_function<double>(erfinv, x); }
237 mat erfinv(const mat &x) { return apply_function<double>(erfinv, x); }
238 
239 // Complementary error function
240 vec erfc(const vec &x) { return apply_function<double>(::erfc, x); }
241 mat erfc(const mat &x) { return apply_function<double>(::erfc, x); }
242 
243 // Q-function
244 vec Qfunc(const vec &x) { return apply_function<double>(Qfunc, x); }
245 mat Qfunc(const mat &x) { return apply_function<double>(Qfunc, x); }
246 
247 } // namespace itpp
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:21 for IT++ by Doxygen 1.8.2