IT++ Logo
error.cpp
Go to the documentation of this file.
1
31#include <itpp/base/itcompat.h>
32
33
34namespace itpp
35{
36
52std::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
88std::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
100std::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
124std::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 */
154std::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
169double 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
217double Qfunc(double x)
218{
219 return (0.5 * ::erfc(x / 1.41421356237310));
220}
221
222
223// Error function
224vec erf(const vec &x) { return apply_function<double>(::erf, x); }
225mat erf(const mat &x) { return apply_function<double>(::erf, x); }
226cvec erf(const cvec &x)
227{
228 return apply_function<std::complex<double> >(erf, x);
229}
230cmat erf(const cmat &x)
231{
232 return apply_function<std::complex<double> >(erf, x);
233}
234
235// Inverse of error function
236vec erfinv(const vec &x) { return apply_function<double>(erfinv, x); }
237mat erfinv(const mat &x) { return apply_function<double>(erfinv, x); }
238
239// Complementary error function
240vec erfc(const vec &x) { return apply_function<double>(::erfc, x); }
241mat erfc(const mat &x) { return apply_function<double>(::erfc, x); }
242
243// Q-function
244vec Qfunc(const vec &x) { return apply_function<double>(Qfunc, x); }
245mat Qfunc(const mat &x) { return apply_function<double>(Qfunc, x); }
246
247} // namespace itpp
Elementary mathematical functions - header file.
Error functions - header file.
double erfinv(double P)
Inverse of error function.
Definition: error.cpp:169
double Qfunc(double x)
Q-function.
Definition: error.cpp:217
std::complex< double > erf(const std::complex< double > &z)
Error function for complex argument.
Definition: error.cpp:154
vec erfc(const vec &x)
Complementary error function.
Definition: error.cpp:240
#define it_error_if(t, s)
Abort if t is true.
Definition: itassert.h:117
vec log(const vec &x)
The natural logarithm of the elements.
Definition: log_exp.h:241
vec exp(const vec &x)
Exp of the elements of a vector x.
Definition: log_exp.h:155
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
double sign(double x)
Signum function.
Definition: elem_math.h:81
IT++ compatibility types and functions.
itpp namespace
Definition: itmex.h:37
std::complex< double > cerfc_continued_fraction(const std::complex< double > &z)
Definition: error.cpp:52
const double pi
Constant Pi.
Definition: misc.h:103
std::complex< double > cerf_rybicki(const std::complex< double > &z)
Definition: error.cpp:124
const double eps
Constant eps.
Definition: misc.h:109
std::complex< double > cerf_series(const std::complex< double > &z)
Definition: error.cpp:100
std::complex< double > cerf_continued_fraction(const std::complex< double > &z)
Complementary function to cerfc_continued_fraction.
Definition: error.cpp:88
bin abs(const bin &inbin)
absolute value of bin
Definition: binary.h:174
int abs(const itpp::bin &inbin)
absolute value of bin
Definition: binary.h:186
SourceForge Logo

Generated on Tue Jan 24 2023 00:00:00 for IT++ by Doxygen 1.9.6