errorfunctions/
complexerrorfunctions.rs

1
2use num::complex::Complex;
3
4use std::f64::consts::PI;
5use crate::realerrorfunctions::RealErrorFunctions;
6use crate::auxilliary::*;
7
8type Complex64 = Complex<f64>;
9
10
11// The Faddeeva function
12
13pub fn w_with_relerror(z: Complex64, mut relerr: f64) -> Complex64 {
14    if z.re == 0.0 {
15        return Complex64::new((z.im).erfcx(), z.re); // give correct sign of 0 in w.im
16    } else if z.im == 0.0 {
17        return Complex64::new((-z.re * z.re).exp(), (z.re).w_im());
18    }
19
20    let a: f64;
21    let a2: f64;
22    let c: f64;
23    if relerr <= f64::EPSILON {
24        relerr = f64::EPSILON;
25        a = 0.518321480430085929872;  // pi / sqrt(-log(eps*0.5))
26        c = 0.329973702884629072537;  // (2/pi) * a;
27        a2 = 0.268657157075235951582; // a^2
28    } else {
29        if relerr > 0.1 {
30            relerr = 0.1; // not sensible to compute < 1 digit
31        }
32        a = PI / (-(relerr * 0.5).log2()).sqrt();
33        c = (2.0 / PI) * a;
34        a2 = a * a;
35    }
36
37    let x: f64 = z.re.abs();
38    let y: f64 = z.im;
39    let ya: f64 = y.abs();
40
41    let ret: Complex64; // return value
42
43    let mut sum1 = 0.0;
44    let mut sum2 = 0.0;
45    let mut sum3 = 0.0;
46    let mut sum4 = 0.0;
47    let mut sum5 = 0.0;
48
49    // Use a continued fraction for large |z|, as it is faster.
50    // As pointed out by M. Zaghloul, the continued fraction seems to give a large relative error in
51    //    Re w(z) for |x| ~ 6 and small |y|, so use algorithm 816 in this region.
52
53    if ya > 7.0 || (x > 6.0 && (ya > 0.1 || (x > 8.0 && ya > 1e-10) || x > 28.0)) {
54        // Poppe & Wijers suggest using a number of terms
55        //      nu = 3 + 1442 / (26*rho + 77)
56        // where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
57        // They only use this expansion for rho >= 1, but rho a little less than 1 seems okay too.
58        // Instead, I did my own fit to a slightly different function that avoids the hypotenuse calculation, using NLopt to minimize
59        // the sum of the squares of the errors in nu with the constraint that the estimated nu be >= minimum nu to attain machine precision.
60        // I also separate the regions where nu == 2 and nu == 1.
61
62        const ISPI: f64 = 0.56418958354775628694807945156; // 1 / sqrt(pi)
63        let xs: f64 = if y < 0.0 { -z.re } else { z.re }; // Compute for -z if y < 0
64        if x + ya > 4000.0 {
65            // nu <= 2
66            if x + ya > 1.0e7 {
67                // nu == 1, w(z) = i/sqrt(pi) / z
68                // scale to avoid overflow
69                if x > ya {
70                    let yax = ya / xs;
71                    let denom = ISPI / (xs + yax * ya);
72                    ret = Complex64::new(denom * yax, denom);
73                } else {
74                    if ya.is_infinite() && ya.is_sign_positive() {
75                        if x.is_nan() || y < 0.0 {
76                            return Complex64::new(f64::NAN, f64::NAN);
77                        } else {
78                            return Complex64::new(0.0, 0.0);
79                        }
80                    } else {
81                        let xya = xs / ya;
82                        let denom = ISPI / (xya * xs + ya);
83                        ret = Complex64::new(denom, denom * xya);
84                    }
85                }
86            } else {
87                // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
88                let dr = xs * xs - ya * ya - 0.5;
89                let di = 2.0 * xs * ya;
90                let denom = ISPI / (dr * dr + di * di);
91                ret = Complex64::new(denom * (xs * di - ya * dr), denom * (xs * dr + ya * di));
92            }
93        } else {
94            // compute nu(z) estimate and do general continued fraction
95            const C0: f64 = 3.9;
96            const C1: f64 = 11.398;
97            const C2: f64 = 0.08254;
98            const C3: f64 = 0.1421;
99            const C4: f64 = 0.2023; // fit
100            let nu = (C0 + C1 / (C2 * x + C3 * ya + C4)).floor();
101            let mut wr = xs;
102            let mut wi = ya;
103            let mut nu = 0.5 * (nu - 1.0);
104            while nu > 0.4 {
105                // w <- z - nu/w:
106                let denom = nu / (wr * wr + wi * wi);
107                wr = xs - wr * denom;
108                wi = ya + wi * denom;
109                nu -= 0.5;
110            }
111            {
112                // w(z) = i/sqrt(pi) / w:
113                let denom = ISPI / (wr * wr + wi * wi);
114                ret = Complex64::new(denom * wi, denom * wr);
115            }
116        }
117
118        if y < 0.0 {
119            // use w(z) = 2.0*exp(-z*z) - w(-z), but be careful of overflow in exp(-z*z) = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
120            return 2.0 * Complex64::new((ya - xs) * (xs + ya), 2.0 * xs * y).exp() - ret;
121        } else {
122            return ret;
123        }
124    } else {
125        // Note: The test that seems to be suggested in the paper is x < sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2)
126        // underflows to zero and sum1,sum2,sum4 are zero.  However, long before this occurs, the sum1,sum2,sum4 contributions are
127        // negligible in double precision; I find that this happens for x > about 6, for all y.  On the other hand, I find that the case
128        // where we compute all of the sums is faster (at least with the precomputed expa2n2 table) until about x=10.  Furthermore, if we
129        // try to compute all of the sums for x > 20, I find that we sometimes run into numerical problems because underflow/overflow
130        // problems start to appear in the various coefficients of the sums, below.  Therefore, we use x < 10 here.
131
132        if x < 10.0 {
133            let mut prod2ax = 1.0;
134            let mut prodm2ax = 1.0;
135            let expx2;
136
137            if y.is_nan() {
138                return Complex64::new(y, y);
139            }
140
141            // Somewhat ugly copy-and-paste duplication here, but I see significant speedups from using the special-case code
142            // with the precomputed exponential, and the x < 5e-4 special case is needed for accuracy.
143
144            if relerr == f64::EPSILON {
145                // Use precomputed exp(-a2*(n*n)) table
146                if x < 5e-4 {
147                    // compute sum4 and sum5 together as sum5-sum4
148                    let x2 = x * x;
149                    expx2 = 1.0 - x2 * (1.0 - 0.5 * x2); // exp(-x*x) via Taylor
150
151                    // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
152                    let ax2 = 1.036642960860171859744 * x; // 2*a*x
153                    let exp2ax =
154                        1.0 + ax2 * (1.0 + ax2 * (0.5 + 0.166666666666666666667 * ax2));
155                    let expm2ax =
156                        1.0 - ax2 * (1.0 - ax2 * (0.5 - 0.166666666666666666667 * ax2));
157                    let mut n: i32 = 1;
158                    loop {
159                        let coef = EXPA2N2[n as usize - 1] * expx2
160                            / (a2 * (f64::from(n) * f64::from(n)) + y * y);
161                        prod2ax *= exp2ax;
162                        prodm2ax *= expm2ax;
163                        sum1 += coef;
164                        sum2 += coef * prodm2ax;
165                        sum3 += coef * prod2ax;
166
167                        // really = sum5 - sum4
168                        sum5 += coef
169                            * (2.0 * a)
170                            * f64::from(n)
171                            * sinh_taylor((2.0 * a) * f64::from(n) * x);
172
173                        // test convergence via sum3
174                        if coef * prod2ax < relerr * sum3 {
175                            break;
176                        }
177                        n += 1;
178                    }
179                } else {
180                    // x > 5e-4, compute sum4 and sum5 separately
181                    expx2 = (-x * x).exp();
182                    let exp2ax = ((2.0 * a) * x).exp();
183                    let expm2ax = 1.0 / exp2ax;
184                    let mut n: i32 = 1;
185                    loop {
186                        let coef = EXPA2N2[n as usize - 1] * expx2
187                            / (a2 * (f64::from(n) * f64::from(n)) + y * y);
188                        prod2ax *= exp2ax;
189                        prodm2ax *= expm2ax;
190                        sum1 += coef;
191                        sum2 += coef * prodm2ax;
192                        sum4 += (coef * prodm2ax) * (a * n as f64);
193                        sum3 += coef * prod2ax;
194                        sum5 += (coef * prod2ax) * (a * n as f64);
195                        // test convergence via sum5, since this sum has the slowest decay
196                        if coef * prod2ax * a * f64::from(n) < relerr * sum5 {
197                            break;
198                        }
199                        n += 1;
200                    }
201                }
202            } else {
203                // relerr != f64::EPSILON, compute exp(-a2*(n*n)) on the fly
204                let exp2ax = (2.0 * a * x).exp();
205                let expm2ax = 1.0 / exp2ax;
206                if x < 5e-4 {
207                    // compute sum4 and sum5 together as sum5-sum4
208                    let x2 = x * x;
209                    expx2 = 1.0 - x2 * (1.0 - 0.5 * x2); // exp(-x*x) via Taylor
210                    let mut n: i32 = 1;
211                    loop {
212                        let coef = (-a2 * f64::from(n) * f64::from(n)).exp() * expx2
213                            / (a2 * (f64::from(n) * f64::from(n)) + y * y);
214                        prod2ax *= exp2ax;
215                        prodm2ax *= expm2ax;
216                        sum1 += coef;
217                        sum2 += coef * prodm2ax;
218                        sum3 += coef * prod2ax;
219
220                        // Really = sum5 - sum4
221                        sum5 += coef
222                            * (2.0 * a)
223                            * f64::from(n)
224                            * sinh_taylor((2.0 * a) * f64::from(n) * x);
225
226                        // Test convergence via sum3
227                        if coef * prod2ax < relerr * sum3 {
228                            break;
229                        }
230                        n += 1;
231                    }
232                } else {
233                    // x > 5e-4, compute sum4 and sum5 separately
234                    expx2 = (-x * x).exp();
235                    let mut n: i32 = 1;
236                    loop {
237                        let coef = (-a2 * (f64::from(n) * f64::from(n))).exp() * expx2
238                            / (a2 * (f64::from(n) * f64::from(n)) + y * y);
239                        prod2ax *= exp2ax;
240                        prodm2ax *= expm2ax;
241                        sum1 += coef;
242                        sum2 += coef * prodm2ax;
243                        sum4 += (coef * prodm2ax) * (a * n as f64);
244                        sum3 += coef * prod2ax;
245                        sum5 += (coef * prod2ax) * (a * n as f64);
246                        // Test convergence via sum5, since this sum has the slowest decay
247                        if (coef * prod2ax) * (a * n as f64) < relerr * sum5 {
248                            break;
249                        }
250                        n += 1;
251                    }
252                }
253            }
254
255            let expx2erfcxy: f64 = if y > -6.0 {
256                expx2 * y.erfcx()
257            } else {
258                2.0 * (y * y - x * x).exp()
259            }; // avoid spurious overflow for large negative y
260            if y > 5.0 {
261                // imaginary terms cancel
262                let sinxy = (x * y).sin();
263                ret = Complex64::new(
264                    (expx2erfcxy - c * y * sum1) * (2.0 * x * y).cos()
265                        + (c * x * expx2) * sinxy * sinc(x * y, sinxy),
266                    0.0,
267                );
268            } else {
269                let xs = z.re;
270                let sinxy = (xs * y).sin();
271                let sin2xy = (2.0 * xs * y).sin();
272                let cos2xy = (2.0 * xs * y).cos();
273                let coef1 = expx2erfcxy - c * y * sum1;
274                let coef2 = c * xs * expx2;
275                ret = Complex64::new(
276                    coef1 * cos2xy + coef2 * sinxy * sinc(xs * y, sinxy),
277                    coef2 * sinc(2.0 * xs * y, sin2xy) - coef1 * sin2xy,
278                );
279            }
280        } else {
281            // x >= 10: only sum3 & sum5 contribute (see above note)
282            if x.is_nan() {
283                return Complex64::new(x, x);
284            }
285            if y.is_nan() {
286                return Complex64::new(y, y);
287            }
288
289            ret = Complex64::new((-x * x).exp(), 0.0); // |y| < 1e-10, so we only need exp(-x*x) term
290
291            // Round instead of ceil as in original paper; note that x/a > 1 here
292            let n0 = (x / a + 0.5).floor(); // sum in both directions, starting at n0
293            let dx = a * n0 - x;
294            sum3 = (-dx * dx).exp() / (a2 * (n0 * n0) + y * y);
295            sum5 = a * n0 * sum3;
296            let exp1 = (4.0 * a * dx).exp();
297            let mut exp1dn = 1.0;
298            let mut dn: i32 = 1;
299            while dn < n0 as i32 {
300                // loop over n0-dn and n0+dn terms
301                let np = n0 + dn as f64;
302                let nm = n0 - dn as f64;
303                let mut tp = (-sqr(a * f64::from(dn) + dx)).exp();
304                exp1dn *= exp1;
305                let mut tm = tp * exp1dn; // trick to get tm from tp
306                tp /= a2 * (np * np) + y * y;
307                tm /= a2 * (nm * nm) + y * y;
308                sum3 += tp + tm;
309                sum5 += a * (np * tp + nm * tm);
310                if a * (np * tp + nm * tm) < relerr * sum5 {
311                    return ret + Complex64::new( (0.5 * c) * y * (sum2 + sum3), (0.5 * c) * (sum5 - sum4).copysign(z.re) );
312                }
313                dn += 1;
314            }
315
316            loop {
317                // loop over n0+dn terms only (since n0-dn <= 0)
318                let np = n0 + dn as f64;
319                dn += 1;
320                let tp = (-sqr(a * f64::from(dn) + dx)).exp() / (a2 * (np * np) + y * y);
321                sum3 += tp;
322                sum5 += a * np * tp;
323                if a * np * tp < relerr * sum5 {
324                    return ret + Complex64::new( (0.5 * c) * y * (sum2 + sum3), (0.5 * c) * (sum5 - sum4).copysign(z.re) );
325                }
326            }
327        }
328    }
329
330    return ret + Complex64::new( (0.5 * c) * y * (sum2 + sum3), (0.5 * c) * (sum5 - sum4).copysign(z.re) );
331}
332
333
334
335
336
337
338
339
340// The complex scaled complementary error function
341
342pub fn erfcx_with_relerror(z: Complex64, relerr: f64) -> Complex64 {
343    w_with_relerror(Complex64::new(-z.im, z.re), relerr)
344}
345
346
347
348
349
350
351
352// The complex error function
353
354pub fn erf_with_relerror(z: Complex64, relerr: f64) -> Complex64 {
355    let x = z.re;
356    let y = z.im;
357
358    // If the imaginary part is 0, make sure the sign of 0 is preserved
359    if y == 0.0 {
360        return Complex64::new(x.erf(), y);
361    }
362
363    // If the real part is 0, make sure the sign of 0 is preserved.
364    // Handle the case of y -> Inf limit manually, since exp(y^2) -> Inf
365    // but Im[w(y)] -> 0, so IEEE will give us a NaN when it should be Inf.
366    if x == 0.0 {
367        let imag = if y * y > 720.0 {
368            if y > 0.0 {
369                f64::INFINITY
370            } else {
371                f64::NEG_INFINITY
372            }
373        } else {
374            (y * y).exp() * y.w_im()
375        };
376
377        return Complex64::new(x, imag);
378    }
379
380    let m_re_z2: f64 = (y - x) * (x + y);                    // Re(-z^2), being careful of overflow
381    let m_im_z2: f64 = -2.0 * x * y;                         // Im(-z^2)
382    if m_re_z2 < -750.0 {
383        // underflow
384        if x >= 0.0 {
385            return Complex64::new(1.0, 0.0);
386        } else {
387            return Complex64::new(-1.0, 0.0);
388        }
389    }
390
391    // Handle positive and negative x via different formulas, using the mirror symmetries of w,
392    // to avoid overflow/underflow problems from multiplying exponentially large and small quantities.
393    if x >= 0.0 {
394        if x < 8.0e-2 {
395            if y.abs() < 1.0e-2 {
396                // Use Taylor series for small |z|, to avoid cancellation inaccuracy
397                //   erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
398
399                let mz2 = Complex64::new(m_re_z2, m_im_z2); // -z^2
400                return z * (1.1283791670955125739 + mz2 * (0.37612638903183752464 + mz2 * (0.11283791670955125739 + mz2 * (0.026866170645131251760 + mz2 * 0.0052239776254421878422))));
401            } else if m_im_z2.abs() < 5.0e-3 && x < 5.0e-3 {
402                // For small |x| and small |xy|, use Taylor series to avoid cancellation inaccuracy:
403                //     erf(x+iy) = erf(iy)
404                //        + 2*exp(y^2)/sqrt(pi) *
405                //          [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
406                //            - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
407                //   where:
408                //      erf(iy) = exp(y^2) * Im[w(y)]
409
410                let x2: f64 = x * x;
411                let y2: f64 = y * y;
412                let expy2: f64 = y2.exp();
413                return Complex64::new(
414                    expy2 * x * (1.1283791670955125739
415                            - x2 * (0.37612638903183752464 + 0.75225277806367504925 * y2)
416                            + x2 * x2 * (0.11283791670955125739 + y2 * (0.45135166683820502956 + 0.15045055561273500986 * y2))),
417                    expy2 * (y.w_im() - x2 * y * (1.1283791670955125739 - x2 * (0.56418958354775628695 + 0.37612638903183752464 * y2)))
418                );
419            }
420        }
421
422        // Don't use the complex exp function, since that will produce spurious NaN values when multiplying w in an overflow situation.
423        return 1.0 - m_re_z2.exp() * Complex64::new(m_im_z2.cos(), m_im_z2.sin()) * w_with_relerror(Complex64::new(-y, x), relerr);
424    } else {
425        // x < 0
426        if x > -8.0e-2 {
427            // duplicate from above to avoid abs(x) call
428            if y.abs() < 1.0e-2 {
429                // Use Taylor series for small |z|, to avoid cancellation inaccuracy
430                //   erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
431
432                let mz2: Complex64 = Complex64::new(m_re_z2, m_im_z2); // -z^2
433                return z * (1.1283791670955125739 + mz2 * (0.37612638903183752464 + mz2 * (0.11283791670955125739 + mz2 * (0.026866170645131251760 + mz2 * 0.0052239776254421878422))));
434            } else if m_im_z2.abs() < 5.0e-3 && x > -5.0e-3 {
435                // For small |x| and small |xy|, use Taylor series to avoid cancellation inaccuracy:
436                //     erf(x+iy) = erf(iy)
437                //        + 2*exp(y^2)/sqrt(pi) *
438                //          [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
439                //            - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
440                //   where:
441                //      erf(iy) = exp(y^2) * Im[w(y)]
442
443                let x2: f64 = x * x;
444                let y2: f64 = y * y;
445                let expy2: f64 = y2.exp();
446                return Complex64::new(
447                    expy2 * x * (1.1283791670955125739
448                            - x2 * (0.37612638903183752464 + 0.75225277806367504925 * y2)
449                            + x2 * x2 * (0.11283791670955125739 + y2 * (0.45135166683820502956 + 0.15045055561273500986 * y2))),
450                    expy2 * (y.w_im() - x2 * y * (1.1283791670955125739 - x2 * (0.56418958354775628695 + 0.37612638903183752464 * y2)))
451                );
452            }
453        } else if x.is_nan() {
454            if y == 0.0 {
455                return Complex64::new(f64::NAN, 0.0);
456            } else {
457                return Complex64::new(f64::NAN, f64::NAN);
458            }
459        }
460
461        // Don't use the complex exp function, since that will produce spurious NaN values when multiplying w in an overflow situation.
462        // Although in principle complex multiplication should be associative, in numerical applications this is not necessarily the case.
463        // In the following I therefore deliberately write a * (b * c) rather than a * b * c. If not, one of the test cases will break
464        // and return (NaN, -Inf) rather than (-Inf, -Inf).
465
466        m_re_z2.exp() * (Complex64::new(m_im_z2.cos(), m_im_z2.sin()) * w_with_relerror(Complex64::new(y, -x), relerr)) - 1.0
467    }
468}
469
470
471
472
473
474
475// The imaginary error function
476
477pub fn erfi_with_relerror(z: Complex64, relerr: f64) -> Complex64 {
478    let e: Complex64 = erf_with_relerror(Complex64::new(-z.im, z.re), relerr);
479    Complex64::new(e.im, -e.re)
480}
481
482
483
484
485
486// The complex complementary error function
487
488pub fn erfc_with_relerror(z: Complex64, relerr: f64) -> Complex64 {
489    let x: f64 = z.re;
490    let y: f64 = z.im;
491
492    if x == 0.0 {
493        // Handle y -> Inf limit manualyy, since exp(y^2) -> Inf but Im[w(y)] -> 0, so
494        // IEEE will give us a NaN when it should be Inf
495
496        if y * y > 720.0 {
497            if y > 0.0 {
498                return Complex64::new(1.0, f64::NEG_INFINITY);
499            } else {
500                return Complex64::new(1.0, f64::INFINITY);
501            }
502        } else {
503            return Complex64::new(1.0, -(y * y).exp() * y.w_im());
504        }
505    }
506
507    if y == 0.0 {
508        if x * x > 750.0 {
509            // underflow
510            if x >= 0.0 {
511                return Complex64::new(0.0, -y); // Preserve sign of 0
512            } else {
513                return Complex64::new(2.0, -y); // Preserve sign of 0
514            }
515        } else {
516            if x >= 0.0 {
517                return Complex64::new((-x*x).exp() * x.erfcx(), -y);
518            } else {
519                return Complex64::new(2.0 - (-x*x).exp() * (-x).erfcx(), -y);
520            }
521        }
522    }
523
524    let m_re_z2: f64 = (y - x) * (x + y);      // Re(-z^2), being careful of overflow
525    let m_im_z2: f64 = -2.0 * x * y;           // Im(-z^2)
526    if m_re_z2 < -750.0 {
527        // underflow
528        if x >= 0.0 {
529            return Complex64::new(0.0, 0.0);
530        } else {
531            return Complex64::new(2.0, 0.0);
532        }
533    }
534
535    if x >= 0.0 {
536        return Complex64::new(m_re_z2, m_im_z2).exp() * w_with_relerror(Complex64::new(-y, x), relerr);
537    } else {
538        return 2.0 - Complex64::new(m_re_z2, m_im_z2).exp() * w_with_relerror(Complex64::new(y, -x), relerr);
539    }
540}
541
542
543
544
545
546// Dawson's function 
547
548pub fn dawson_with_relerror(z: Complex64, relerr: f64) -> Complex64 {
549    const SPI2: f64 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
550    let x = z.re;
551    let y = z.im;
552
553    // Handle axes separately for speed & proper handling of x or y = Inf or NaN
554
555    if y == 0.0 {
556        return Complex64::new(SPI2 * x.w_im(), -y); // preserve sign of 0
557    }
558
559    if x == 0.0 {
560        let y2 = y * y;
561        if y2 < 2.5e-5 {
562            // Taylor expansion
563            return Complex64::new( x, y * (1. + y2 * (0.6666666666666666666666666666666666666667 + y2 * 0.26666666666666666666666666666666666667)));
564        } else {
565            // Preserve sign of 0
566            if y >= 0.0 {
567                return Complex64::new(x, SPI2 * (y2.exp() - y.erfcx()));
568            } else {
569                return Complex64::new(x, SPI2 * ((-y).erfcx() - y2.exp()));
570            }
571        }
572    }
573
574    let m_re_z2 = (y - x) * (x + y);             // Re(-z^2), being careful of overflow
575    let m_im_z2 = -2.0 * x * y;                  // Im(-z^2)
576    let mz2 = Complex64::new(m_re_z2, m_im_z2);  // -z^2
577
578    // Handle positive and negative x via different formulas, using the mirror symmetries of w, to avoid overflow/underflow
579    // problems from multiplying exponentially large and small quantities.
580
581    if y >= 0.0 {
582        if y < 5e-3 {
583            if x.abs() < 5e-3 {
584                // Use Taylor series for small |z|, to avoid cancellation inaccuracy: dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
585                return z * (1. + mz2 * (0.6666666666666666666666666666666666666667 + mz2 * 0.2666666666666666666666666666666666666667));
586            } else if m_im_z2.abs() < 5.0e-3 {
587                // (**) For small |y| and small |xy|, use Taylor series to avoid cancellation inaccuracy:
588                //     dawson(x + iy) = D + y^2 (D + x - 2Dx^2) + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
589                //                        + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
590                //                        + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
591                // where D = dawson(x)
592                //
593                // However, for large |x|, 2Dx -> 1 which gives cancellation problems in this series (many of the leading terms cancel).
594                // So, for large |x|, we need to substitute a continued-fraction expansion for D:
595                //     dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))
596                //
597                // The 6 terms shown here seems to be the minimum needed to be accurate as soon as the simpler Taylor expansion above starts
598                // breaking down.  Using this 6-term expansion, factoring out the denominator, and simplifying with Maple, we obtain:
599                //     Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
600                //     Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4
601                //
602                // Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction expansion for the real part,
603                // and a 2-term expansion for the imaginary part. This avoids overflow problems for huge |x|.  This yields:
604                //     Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
605                //     Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)
606
607                let x2 = x * x;
608                if x2 > 1600.0 {
609                    // |x| > 40
610                    let y2 = y * y;
611                    if x2 > 25.0e14 {
612                        // |x| > 5e7
613                        let xy2 = (x * y) * (x * y);
614                        return Complex64::new(
615                            (0.5 + y2 * (0.5 + 0.25 * y2 - 0.16666666666666666667 * xy2)) / x,
616                            y * (-1.0
617                                + y2 * (-0.66666666666666666667
618                                    + 0.13333333333333333333 * xy2
619                                    - 0.26666666666666666667 * y2))
620                                / (2.0 * x2 - 1.0),
621                        );
622                    }
623                    return (1. / (-15. + x2 * (90. + x2 * (-60. + 8. * x2))))
624                        * Complex64::new(
625                            x * (33. + x2 * (-28. + 4. * x2) + y2 * (18. - 4. * x2 + 4. * y2)),
626                            y * (-15. + x2 * (24. - 4. * x2) + y2 * (4. * x2 - 10. - 4. * y2)),
627                        );
628                } else {
629                    let d = SPI2 * x.w_im();
630                    let y2 = y * y;
631                    return Complex64::new(
632                        d + y2 * (d + x - 2. * d * x2)
633                            + y2 * y2
634                                * (d * (0.5 - x2 * (2. - 0.66666666666666666667 * x2))
635                                    + x * (0.83333333333333333333
636                                        - 0.33333333333333333333 * x2)),
637                        y * (1. - 2. * d * x
638                            + y2 * 0.66666666666666666667 * (1. - x2 - d * x * (3. - 2. * x2))
639                            + y2 * y2
640                                * (0.26666666666666666667
641                                    - x2 * (0.6 - 0.13333333333333333333 * x2)
642                                    - d * x
643                                        * (1.
644                                            - x2 * (1.3333333333333333333
645                                                - 0.26666666666666666667 * x2)))),
646                    );
647                }
648            }
649        }
650        let res: Complex64 = mz2.exp() - w_with_relerror(z, relerr);
651        return SPI2 * Complex64::new(-res.im, res.re);
652    } else {
653        // y < 0
654        if y > -5.0e-3 {
655            // duplicate from above to avoid abs(x) call
656            if x.abs() < 5.0e-3 {
657                // Use Taylor series for small |z|, to avoid cancellation inaccuracy: dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
658                return z * (1. + mz2 * (0.6666666666666666666666666666666666666667 + mz2 * 0.2666666666666666666666666666666666666667));
659            } else if m_im_z2.abs() < 5.0e-3 {
660                // For small |y| and small |xy|, use Taylor series to avoid cancellation inaccuracy.
661                // See explanation above at (**).
662                let x2 = x * x;
663                if x2 > 1600.0 {
664                    // |x| > 40
665                    let y2 = y * y;
666                    if x2 > 25.0e14 {
667                        // |x| > 5e7
668                        let xy2 = (x * y) * (x * y);
669                        return Complex64::new(
670                            (0.5 + y2 * (0.5 + 0.25 * y2 - 0.16666666666666666667 * xy2)) / x,
671                            y * (-1.0
672                                + y2 * (-0.66666666666666666667
673                                    + 0.13333333333333333333 * xy2
674                                    - 0.26666666666666666667 * y2))
675                                / (2.0 * x2 - 1.0),
676                        );
677                    }
678                    return (1. / (-15. + x2 * (90. + x2 * (-60. + 8. * x2))))
679                        * Complex64::new(
680                            x * (33. + x2 * (-28. + 4. * x2) + y2 * (18. - 4. * x2 + 4. * y2)),
681                            y * (-15. + x2 * (24. - 4. * x2) + y2 * (4. * x2 - 10. - 4. * y2))
682                        );
683                } else {
684                    let d = SPI2 * x.w_im();
685                    let y2 = y * y;
686                    return Complex64::new(
687                        d + y2 * (d + x - 2. * d * x2)
688                            + y2 * y2
689                                * (d * (0.5 - x2 * (2. - 0.66666666666666666667 * x2))
690                                    + x * (0.83333333333333333333
691                                        - 0.33333333333333333333 * x2)),
692                        y * (1. - 2. * d * x
693                            + y2 * 0.66666666666666666667 * (1. - x2 - d * x * (3. - 2. * x2))
694                            + y2 * y2
695                                * (0.26666666666666666667
696                                    - x2 * (0.6 - 0.13333333333333333333 * x2)
697                                    - d * x
698                                        * (1.
699                                            - x2 * (1.3333333333333333333
700                                                - 0.26666666666666666667 * x2)))),
701                    );
702                }
703            }
704        } else if y.is_nan() {
705            if x == 0.0 {
706                return Complex64::new(0.0, f64::NAN);
707            } else {
708                return Complex64::new(f64::NAN, f64::NAN);
709            }
710        }
711        let res: Complex64 = w_with_relerror(-z, relerr) - mz2.exp();
712        return SPI2 * Complex64::new(-res.im, res.re);
713    }
714}
715
716
717
718
719
720pub trait ComplexErrorFunctions {
721    fn erfcx(self) -> Self;
722    fn erf(self) -> Self;
723    fn w(self) -> Self;
724    fn erfi(self) -> Self;      
725    fn erfc(self) -> Self;
726    fn dawson(self) -> Self;
727}
728
729
730
731impl ComplexErrorFunctions for Complex64 {
732
733    fn erfcx(self) -> Self {
734        erfcx_with_relerror(self, 0.0)
735    }
736
737    fn erf(self) -> Self {
738        erf_with_relerror(self, 0.0)
739    }
740
741    fn w(self) -> Self {
742        w_with_relerror(self, 0.0)
743    }
744
745    fn erfi(self) -> Self {
746        erfi_with_relerror(self, 0.0)
747    }
748
749    fn erfc(self) -> Self {
750        erfc_with_relerror(self, 0.0)
751    }
752
753    fn dawson(self) -> Self {
754        dawson_with_relerror(self, 0.0)
755    }
756}
757
758
759
760
761
762
763
764
765
766#[cfg(test)]
767mod tests {
768
769    use super::*;
770    use crate::auxilliary::relerr;
771
772    #[test]
773    fn test_w_complex() {
774        // Function arguments z.
775
776        let z: [Complex64; 57] = [
777            Complex64::new(624.2, -0.26123),
778            Complex64::new(-0.4, 3.),
779            Complex64::new(0.6, 2.),
780            Complex64::new(-1., 1.),
781            Complex64::new(-1., -9.),
782            Complex64::new(-1., 9.),
783            Complex64::new(-0.0000000234545, 1.1234),
784            Complex64::new(-3., 5.1),
785            Complex64::new(-53.0, 30.1),
786            Complex64::new(0.0, 0.12345),
787            Complex64::new(11.0, 1.0),
788            Complex64::new(-22.0, -2.0),
789            Complex64::new(9.0, -28.0),
790            Complex64::new(21.0, -33.0),
791            Complex64::new(1e5, 1e5),
792            Complex64::new(1e14, 1e14),
793            Complex64::new(-3001.0, -1000.0),
794            Complex64::new(1e160, -1e159),
795            Complex64::new(-6.01, 0.01),
796            Complex64::new(-0.7, -0.7),
797            Complex64::new(2.611780000000000e+01, 4.540909610972489e+03),
798            Complex64::new(0.8e7, 0.3e7),
799            Complex64::new(-20.0, -19.8081),
800            Complex64::new(1e-16, -1.1e-16),
801            Complex64::new(2.3e-8, 1.3e-8),
802            Complex64::new(6.3, -1e-13),
803            Complex64::new(6.3, 1e-20),
804            Complex64::new(1e-20, 6.3),
805            Complex64::new(1e-20, 16.3),
806            Complex64::new(9.0, 1e-300),
807            Complex64::new(6.01, 0.11),
808            Complex64::new(8.01, 1.01e-10),
809            Complex64::new(28.01, 1e-300),
810            Complex64::new(10.01, 1e-200),
811            Complex64::new(10.01, -1e-200),
812            Complex64::new(10.01, 0.99e-10),
813            Complex64::new(10.01, -0.99e-10),
814            Complex64::new(1e-20, 7.01),
815            Complex64::new(-1.0, 7.01),
816            Complex64::new(5.99, 7.01),
817            Complex64::new(1.0, 0.0),
818            Complex64::new(55.0, 0.0),
819            Complex64::new(-0.1, 0.0),
820            Complex64::new(1e-20, 0.0),
821            Complex64::new(0.0, 5e-14),
822            Complex64::new(0.0, 51.0),
823            Complex64::new(f64::INFINITY, 0.0),
824            Complex64::new(f64::NEG_INFINITY, 0.0),
825            Complex64::new(0.0, f64::INFINITY),
826            Complex64::new(0.0, f64::NEG_INFINITY),
827            Complex64::new(f64::INFINITY, f64::INFINITY),
828            Complex64::new(f64::INFINITY, f64::NEG_INFINITY),
829            Complex64::new(f64::NAN, f64::NAN),
830            Complex64::new(f64::NAN, 0.0),
831            Complex64::new(0.0, f64::NAN),
832            Complex64::new(f64::NAN, f64::INFINITY),
833            Complex64::new(f64::INFINITY, f64::NAN),
834        ];
835
836        // w(z) computed with WolframAlpha. WolframAlpha is problematic for some of the above inputs, so I had to use the
837        // continued-fraction expansion in WolframAlpha in some cases, or switch to Maple
838
839        let expected_w: [Complex64; 57] = [
840            Complex64::new( -3.78270245518980507452677445620103199303131110e-7, 0.000903861276433172057331093754199933411710053155),
841            Complex64::new( 0.1764906227004816847297495349730234591778719532788, -0.02146550539468457616788719893991501311573031095617),
842            Complex64::new( 0.2410250715772692146133539023007113781272362309451, 0.06087579663428089745895459735240964093522265589350),
843            Complex64::new( 0.30474420525691259245713884106959496013413834051768, -0.20821893820283162728743734725471561394145872072738),
844            Complex64::new( 7.317131068972378096865595229600561710140617977e34, 8.321873499714402777186848353320412813066170427e34),
845            Complex64::new( 0.0615698507236323685519612934241429530190806818395, -0.00676005783716575013073036218018565206070072304635),
846            Complex64::new( 0.3960793007699874918961319170187598400134746631, -5.593152259116644920546186222529802777409274656e-9),
847            Complex64::new( 0.08217199226739447943295069917990417630675021771804, -0.04701291087643609891018366143118110965272615832184),
848            Complex64::new( 0.00457246000350281640952328010227885008541748668738, -0.00804900791411691821818731763401840373998654987934),
849            Complex64::new(0.8746342859608052666092782112565360755791467973338452, 0.),
850            Complex64::new( 0.00468190164965444174367477874864366058339647648741, 0.0510735563901306197993676329845149741675029197050),
851            Complex64::new( -0.0023193175200187620902125853834909543869428763219, -0.025460054739731556004902057663500272721780776336),
852            Complex64::new( 9.11463368405637174660562096516414499772662584e304, 3.97101807145263333769664875189354358563218932e305),
853            Complex64::new( -4.4927207857715598976165541011143706155432296e281, -2.8019591213423077494444700357168707775769028e281),
854            Complex64::new( 2.820947917809305132678577516325951485807107151e-6, 2.820947917668257736791638444590253942253354058e-6),
855            Complex64::new( 2.82094791773878143474039725787438662716372268e-15, 2.82094791773878143474039725773333923127678361e-15),
856            Complex64::new( -0.0000563851289696244350147899376081488003110150498, -0.000169211755126812174631861529808288295454992688),
857            Complex64::new( -5.586035480670854326218608431294778077663867e-162, 5.586035480670854326218608431294778077663867e-161),
858            Complex64::new( 0.00016318325137140451888255634399123461580248456, -0.095232456573009287370728788146686162555021209999),
859            Complex64::new( 0.69504753678406939989115375989939096800793577783885, -1.8916411171103639136680830887017670616339912024317),
860            Complex64::new( 0.0001242418269653279656612334210746733213167234822, 7.145975826320186888508563111992099992116786763e-7),
861            Complex64::new( 2.318587329648353318615800865959225429377529825e-8, 6.182899545728857485721417893323317843200933380e-8),
862            Complex64::new( -0.0133426877243506022053521927604277115767311800303, -0.0148087097143220769493341484176979826888871576145),
863            Complex64::new( 1.00000000000000012412170838050638522857747934, 1.12837916709551279389615890312156495593616433e-16),
864            Complex64::new( 0.9999999853310704677583504063775310832036830015, 2.595272024519678881897196435157270184030360773e-8),
865            Complex64::new( -1.4731421795638279504242963027196663601154624e-15, 0.090727659684127365236479098488823462473074709),
866            Complex64::new( 5.79246077884410284575834156425396800754409308e-18, 0.0907276596841273652364790985059772809093822374),
867            Complex64::new( 0.0884658993528521953466533278764830881245144368, 1.37088352495749125283269718778582613192166760e-22),
868            Complex64::new( 0.0345480845419190424370085249304184266813447878, 2.11161102895179044968099038990446187626075258e-23),
869            Complex64::new( 6.63967719958073440070225527042829242391918213e-36, 0.0630820900592582863713653132559743161572639353),
870            Complex64::new( 0.00179435233208702644891092397579091030658500743634, 0.0951983814805270647939647438459699953990788064762),
871            Complex64::new( 9.09760377102097999924241322094863528771095448e-13, 0.0709979210725138550986782242355007611074966717),
872            Complex64::new( 7.2049510279742166460047102593255688682910274423e-304, 0.0201552956479526953866611812593266285000876784321),
873            Complex64::new( 3.04543604652250734193622967873276113872279682e-44, 0.0566481651760675042930042117726713294607499165),
874            Complex64::new( 3.04543604652250734193622967873276113872279682e-44, 0.0566481651760675042930042117726713294607499165),
875            Complex64::new( 0.5659928732065273429286988428080855057102069081e-12, 0.056648165176067504292998527162143030538756683302),
876            Complex64::new( -0.56599287320652734292869884280802459698927645e-12, 0.0566481651760675042929985271621430305387566833029),
877            Complex64::new( 0.0796884251721652215687859778119964009569455462, 1.11474461817561675017794941973556302717225126e-22),
878            Complex64::new( 0.07817195821247357458545539935996687005781943386550, -0.01093913670103576690766705513142246633056714279654),
879            Complex64::new( 0.04670032980990449912809326141164730850466208439937, 0.03944038961933534137558064191650437353429669886545),
880            Complex64::new( 0.36787944117144232159552377016146086744581113103176, 0.60715770584139372911503823580074492116122092866515),
881            Complex64::new(0.0, 0.010259688805536830986089913987516716056946786526145), 
882            Complex64::new( 0.99004983374916805357390597718003655777207908125383, -0.11208866436449538036721343053869621153527769495574),
883            Complex64::new( 0.99999999999999999999999999999999999999990000, 1.12837916709551257389615890312154517168802603e-20),
884            Complex64::new(0.999999999999943581041645226871305192054749891144158, 0.0),
885            Complex64::new(0.0110604154853277201542582159216317923453996211744250, 0.0),
886            Complex64::new(0.0, 0.0),
887            Complex64::new(0.0, 0.0),
888            Complex64::new(0.0, 0.0),
889            Complex64::new(f64::INFINITY, 0.0),
890            Complex64::new(0.0, 0.0),
891            Complex64::new(f64::NAN, f64::NAN),
892            Complex64::new(f64::NAN, f64::NAN),
893            Complex64::new(f64::NAN, f64::NAN),
894            Complex64::new(f64::NAN, 0.0),
895            Complex64::new(f64::NAN, f64::NAN),
896            Complex64::new(f64::NAN, f64::NAN),
897        ];
898
899        let tolerance = 1.0e-13;
900        for n in 0..z.len() {
901            let computed_w = z[n].w();
902            // println!( "{}: z={:?} --> w(z) computed: {:?}  vs expected:  {:?}", n, z[n], computed_w, expected_w[n]);
903            let relative_error_re = relerr(computed_w.re, expected_w[n].re);
904            let relative_error_im = relerr(computed_w.im, expected_w[n].im);
905            assert!(relative_error_re < tolerance);
906            assert!(relative_error_im < tolerance);
907        }
908    }
909
910
911
912    #[test]
913    fn test_erf_complex() {
914
915        let z: [Complex64; 41] = [
916            Complex64::new(1.0, 2.0),
917            Complex64::new(-1.0, 2.0),
918            Complex64::new(1.0, -2.0),
919            Complex64::new(-1.0, -2.0),
920            Complex64::new(9.0, -28.0),
921            Complex64::new(21.0, -33.0),
922            Complex64::new(1.0e3, 1.0e3),
923            Complex64::new(-3001.0, -1000.0),
924            Complex64::new(1.0e160, -1.0e159),
925            Complex64::new(5.1e-3, 1.0e-8),
926            Complex64::new(-4.9e-3, 4.95e-3),
927            Complex64::new(4.9e-3, 0.5),
928            Complex64::new(4.9e-4, -0.5e1),
929            Complex64::new(-4.9e-5, -0.5e2),
930            Complex64::new(5.1e-3, 0.5),
931            Complex64::new(5.1e-4, -0.5e1),
932            Complex64::new(-5.1e-5, -0.5e2),
933            Complex64::new(1e-6, 2e-6),
934            Complex64::new(0.0, 2e-6),
935            Complex64::new(0.0, 2.0),
936            Complex64::new(0.0, 20.0),
937            Complex64::new(0.0, 200.0),
938            Complex64::new(f64::INFINITY, 0.0),
939            Complex64::new(f64::NEG_INFINITY, 0.0),
940            Complex64::new(0.0, f64::INFINITY),
941            Complex64::new(0.0, f64::NEG_INFINITY),
942            Complex64::new(f64::INFINITY, f64::INFINITY),
943            Complex64::new(f64::INFINITY, f64::NEG_INFINITY),
944            Complex64::new(f64::NAN, f64::NAN),
945            Complex64::new(f64::NAN, 0.0),
946            Complex64::new(0.0, f64::NAN),
947            Complex64::new(f64::NAN, f64::INFINITY),
948            Complex64::new(f64::INFINITY, f64::NAN),
949            Complex64::new(1e-3, f64::NAN),
950            Complex64::new(7e-2, 7e-2),
951            Complex64::new(7e-2, -7e-4),
952            Complex64::new(-9e-2, 7e-4),
953            Complex64::new(-9e-2, 9e-2),
954            Complex64::new(-7e-4, 9e-2),
955            Complex64::new(7e-2, 0.9e-2),
956            Complex64::new(7e-2, 1.1e-2)
957        ];
958        
959        let expected_erf: [Complex64; 41] = [
960            Complex64::new(-0.5366435657785650339917955593141927494421, -5.049143703447034669543036958614140565553),
961            Complex64::new(0.5366435657785650339917955593141927494421, -5.049143703447034669543036958614140565553),
962            Complex64::new(-0.5366435657785650339917955593141927494421, 5.049143703447034669543036958614140565553),
963            Complex64::new(0.5366435657785650339917955593141927494421, 5.049143703447034669543036958614140565553),
964            Complex64::new(0.3359473673830576996788000505817956637777e304, -0.1999896139679880888755589794455069208455e304),
965            Complex64::new(0.3584459971462946066523939204836760283645e278, 0.3818954885257184373734213077678011282505e280),
966            Complex64::new(0.9996020422657148639102150147542224526887, 0.00002801044116908227889681753993542916894856),
967            Complex64::new(-1.0, 0.0),
968            Complex64::new(1.0, 0.0),
969            Complex64::new(0.005754683859034800134412990541076554934877, 0.1128349818335058741511924929801267822634e-7),
970            Complex64::new(-0.005529149142341821193633460286828381876955, 0.005585388387864706679609092447916333443570),
971            Complex64::new(0.007099365669981359632319829148438283865814, 0.6149347012854211635026981277569074001219),
972            Complex64::new(0.3981176338702323417718189922039863062440e8, -0.8298176341665249121085423917575122140650e10),
973            Complex64::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
974            Complex64::new(0.007389128308257135427153919483147229573895, 0.6149332524601658796226417164791221815139),
975            Complex64::new(0.4143671923267934479245651547534414976991e8, -0.8298168216818314211557046346850921446950e10),
976            Complex64::new(f64::NEG_INFINITY, f64::NEG_INFINITY),
977            Complex64::new(0.1128379167099649964175513742247082845155e-5, 0.2256758334191777400570377193451519478895e-5),
978            Complex64::new(0.0, 0.2256758334194034158904576117253481476197e-5),
979            Complex64::new(0.0, 18.56480241457555259870429191324101719886),
980            Complex64::new(0.0, 0.1474797539628786202447733153131835124599e173),
981            Complex64::new(0.0, f64::INFINITY),
982            Complex64::new(1.0,0.0),
983            Complex64::new(-1.0,0.0),
984            Complex64::new(0.0,f64::INFINITY),
985            Complex64::new(0.0,f64::NEG_INFINITY),
986            Complex64::new(f64::NAN,f64::NAN),
987            Complex64::new(f64::NAN,f64::NAN),
988            Complex64::new(f64::NAN,f64::NAN),
989            Complex64::new(f64::NAN,0.0),
990            Complex64::new(0.0,f64::NAN),
991            Complex64::new(f64::NAN,f64::NAN),
992            Complex64::new(f64::NAN,f64::NAN),
993            Complex64::new(f64::NAN,f64::NAN),
994            Complex64::new(0.07924380404615782687930591956705225541145, 0.07872776218046681145537914954027729115247),
995            Complex64::new(0.07885775828512276968931773651224684454495, -0.0007860046704118224342390725280161272277506),
996            Complex64::new(-0.1012806432747198859687963080684978759881, 0.0007834934747022035607566216654982820299469),
997            Complex64::new(-0.1020998418798097910247132140051062512527, 0.1010030778892310851309082083238896270340),
998            Complex64::new(-0.0007962891763147907785684591823889484764272, 0.1018289385936278171741809237435404896152),
999            Complex64::new(0.07886408666470478681566329888615410479530, 0.01010604288780868961492224347707949372245),
1000            Complex64::new(0.07886723099940260286824654364807981336591, 0.01235199327873258197931147306290916629654)
1001        ];
1002
1003        let tolerance = 1.0e-13;
1004        for n in 0..z.len() {
1005            let computed_erf = z[n].erf();
1006            let relative_error_re = relerr(computed_erf.re, expected_erf[n].re);
1007            let relative_error_im = relerr(computed_erf.im, expected_erf[n].im);
1008            assert!(relative_error_re < tolerance);
1009            assert!(relative_error_im < tolerance);
1010        }
1011    }
1012
1013
1014
1015    #[test]
1016    fn test_erfi_complex() {
1017        let tolerance = 1.0e-13;
1018        let z = Complex64::new(1.234, 0.5678);
1019        let computed_erfi = erfi_with_relerror(z, 0.0);
1020        let expected_erfi = Complex64::new(1.081032284405373149432716643834106923212, 1.926775520840916645838949402886591180834);
1021        let relative_error_re = relerr(computed_erfi.re, expected_erfi.re);
1022        let relative_error_im = relerr(computed_erfi.im, expected_erfi.im);
1023        assert!(relative_error_re < tolerance);
1024        assert!(relative_error_im < tolerance);
1025    }
1026
1027
1028
1029    #[test]
1030    fn test_erfcx_complex() {
1031        let tolerance = 1.0e-13;
1032        let z = Complex64::new(1.234, 0.5678);
1033        let computed_erfcx = erfcx_with_relerror(z, 0.0);
1034        let expected_erfcx = Complex64::new(0.3382187479799972294747793561190487832579, -0.1116077470811648467464927471872945833154);
1035        let relative_error_re = relerr(computed_erfcx.re, expected_erfcx.re);
1036        let relative_error_im = relerr(computed_erfcx.im, expected_erfcx.im);
1037        assert!(relative_error_re < tolerance);
1038        assert!(relative_error_im < tolerance);
1039    }
1040
1041
1042
1043
1044    #[test]
1045    fn test_erfc_complex() {
1046
1047        let z: [Complex64; 30] = [
1048            Complex64::new(1.0, 2.0),
1049            Complex64::new(-1.0, 2.0),
1050            Complex64::new(1.0, -2.0),
1051            Complex64::new(-1.0, -2.0),
1052            Complex64::new(9.0, -28.0),
1053            Complex64::new(21.0, -33.0),
1054            Complex64::new(1.0e3, 1.0e3),
1055            Complex64::new(-3001.0, -1000.0),
1056            Complex64::new(1.0e160, -1.0e159),
1057            Complex64::new(5.1e-3, 1.0e-8),
1058            Complex64::new(0.0, 2.0e-6),
1059            Complex64::new(0.0, 2.0),
1060            Complex64::new(0.0, 20.0),
1061            Complex64::new(0.0, 200.0),
1062            Complex64::new(2.0e-6, 0.0),
1063            Complex64::new(2.0, 0.0),
1064            Complex64::new(20.0, 0.0),
1065            Complex64::new(200.0, 0.0),
1066            Complex64::new(f64::INFINITY, 0.0),
1067            Complex64::new(f64::NEG_INFINITY, 0.0),
1068            Complex64::new(0.0, f64::INFINITY),
1069            Complex64::new(0.0,f64::NEG_INFINITY),
1070            Complex64::new(f64::INFINITY, f64::INFINITY),
1071            Complex64::new(f64::INFINITY, f64::NEG_INFINITY),
1072            Complex64::new(f64::NAN, f64::NAN),
1073            Complex64::new(f64::NAN, 0.0),
1074            Complex64::new(0.0, f64::NAN),
1075            Complex64::new(f64::NAN, f64::INFINITY),
1076            Complex64::new(f64::INFINITY, f64::NAN),
1077            Complex64::new(88.0, 0.0)
1078        ];
1079        
1080        let expected_erfc: [Complex64; 30] = [
1081            Complex64::new(1.536643565778565033991795559314192749442, 5.049143703447034669543036958614140565553),
1082            Complex64::new(0.4633564342214349660082044406858072505579, 5.049143703447034669543036958614140565553),
1083            Complex64::new(1.536643565778565033991795559314192749442, -5.049143703447034669543036958614140565553),
1084            Complex64::new(0.4633564342214349660082044406858072505579, -5.049143703447034669543036958614140565553),
1085            Complex64::new(-0.3359473673830576996788000505817956637777e304, 0.1999896139679880888755589794455069208455e304),
1086            Complex64::new(-0.3584459971462946066523939204836760283645e278, -0.3818954885257184373734213077678011282505e280),
1087            Complex64::new(0.0003979577342851360897849852457775473112748, -0.00002801044116908227889681753993542916894856),
1088            Complex64::new(2.0, 0.0),
1089            Complex64::new(0.0, 0.0),
1090            Complex64::new(0.9942453161409651998655870094589234450651, -0.1128349818335058741511924929801267822634e-7),
1091            Complex64::new(1.0, -0.2256758334194034158904576117253481476197e-5),
1092            Complex64::new(1.0, -18.56480241457555259870429191324101719886),
1093            Complex64::new(1.0, -0.1474797539628786202447733153131835124599e173),
1094            Complex64::new(1.0, f64::NEG_INFINITY),
1095            Complex64::new(0.9999977432416658119838633199332831406314, 0.0),
1096            Complex64::new(0.004677734981047265837930743632747071389108, 0.0),
1097            Complex64::new(0.5395865611607900928934999167905345604088e-175, 0.0),
1098            Complex64::new(0.0, 0.0),
1099            Complex64::new(0.0, 0.0),
1100            Complex64::new(2.0, 0.0),
1101            Complex64::new(1.0, f64::NEG_INFINITY),
1102            Complex64::new(1.0, f64::INFINITY),
1103            Complex64::new(f64::NAN, f64::NAN),
1104            Complex64::new(f64::NAN, f64::NAN),
1105            Complex64::new(f64::NAN, f64::NAN),
1106            Complex64::new(f64::NAN, 0.0),
1107            Complex64::new(1.0, f64::NAN),
1108            Complex64::new(f64::NAN, f64::NAN),
1109            Complex64::new(f64::NAN, f64::NAN),
1110            Complex64::new(0.0, 0.0)
1111        ];
1112
1113        let tolerance = 1.0e-13;
1114        for n in 0..z.len() {
1115            let computed_erfc = erfc_with_relerror(z[n], 0.);
1116            let relative_error_re = relerr(computed_erfc.re, expected_erfc[n].re);
1117            let relative_error_im = relerr(computed_erfc.im, expected_erfc[n].im);
1118            assert!(relative_error_re < tolerance);
1119            assert!(relative_error_im < tolerance);
1120        }
1121    }
1122
1123
1124
1125
1126    #[test]
1127    fn test_dawson_complex() {
1128        let z: [Complex64; 48] = [
1129            Complex64::new(2.0, 1.0),
1130            Complex64::new(-2.0, 1.0),
1131            Complex64::new(2.0, -1.0),
1132            Complex64::new(-2.0, -1.0),
1133            Complex64::new(-28.0, 9.0),
1134            Complex64::new(33.0, -21.0),
1135            Complex64::new(1.0e3, 1.0e3),
1136            Complex64::new(-1000.0, -3001.0),
1137            Complex64::new(1.0e-8, 5.1e-3),
1138            Complex64::new(4.95e-3, -4.9e-3),
1139            Complex64::new(5.1e-3, 5.1e-3),
1140            Complex64::new(0.5, 4.9e-3),
1141            Complex64::new(-0.5e1, 4.9e-4),
1142            Complex64::new(-0.5e2, -4.9e-5),
1143            Complex64::new(0.5e3, 4.9e-6),
1144            Complex64::new(0.5, 5.1e-3),
1145            Complex64::new(-0.5e1, 5.1e-4),
1146            Complex64::new(-0.5e2, -5.1e-5),
1147            Complex64::new(1e-6, 2e-6),
1148            Complex64::new(2e-6, 0.0),
1149            Complex64::new(2.0, 0.0),
1150            Complex64::new(20.0, 0.0),
1151            Complex64::new(200.0, 0.0),
1152            Complex64::new(0.0, 4.9e-3),
1153            Complex64::new(0.0, -5.1e-3),
1154            Complex64::new(0.0, 2e-6),
1155            Complex64::new(0.0, -2.0),
1156            Complex64::new(0.0, 20.0),
1157            Complex64::new(0.0, -200.0),
1158            Complex64::new(f64::INFINITY, 0.0),
1159            Complex64::new(f64::NEG_INFINITY, 0.0),
1160            Complex64::new(0.0, f64::INFINITY),
1161            Complex64::new(0.0, f64::NEG_INFINITY),
1162            Complex64::new(f64::INFINITY, f64::INFINITY),
1163            Complex64::new(f64::INFINITY,f64::NEG_INFINITY),
1164            Complex64::new(f64::NAN,f64::NAN),
1165            Complex64::new(f64::NAN, 0.0),
1166            Complex64::new(0.0, f64::NAN),
1167            Complex64::new(f64::NAN, f64::INFINITY),
1168            Complex64::new(f64::INFINITY, f64::NAN),
1169            Complex64::new(39.0, 6.4e-5),
1170            Complex64::new(41.0, 6.09e-5),
1171            Complex64::new(4.9e7, 5.0e-11),
1172            Complex64::new(5.1e7, 4.8e-11),
1173            Complex64::new(1.0e9, 2.4e-12),
1174            Complex64::new(1.0e11, 2.4e-14),
1175            Complex64::new(1.0e13, 2.4e-16),
1176            Complex64::new(1.0e300, 2.4e-303)
1177        ];
1178
1179        let expected_dawson: [Complex64; 48] = [
1180            Complex64::new(0.1635394094345355614904345232875688576839, -0.1531245755371229803585918112683241066853),
1181            Complex64::new(-0.1635394094345355614904345232875688576839, -0.1531245755371229803585918112683241066853),
1182            Complex64::new(0.1635394094345355614904345232875688576839, 0.1531245755371229803585918112683241066853),
1183            Complex64::new(-0.1635394094345355614904345232875688576839, 0.1531245755371229803585918112683241066853),
1184            Complex64::new(-0.01619082256681596362895875232699626384420, -0.005210224203359059109181555401330902819419),
1185            Complex64::new(0.01078377080978103125464543240346760257008, 0.006866888783433775382193630944275682670599),
1186            Complex64::new(-0.5808616819196736225612296471081337245459, 0.6688593905505562263387760667171706325749),
1187            Complex64::new(f64::INFINITY, f64::NEG_INFINITY),
1188            Complex64::new(0.1000052020902036118082966385855563526705e-7, 0.005100088434920073153418834680320146441685),
1189            Complex64::new(0.004950156837581592745389973960217444687524, -0.004899838305155226382584756154100963570500),
1190            Complex64::new(0.005100176864319675957314822982399286703798, 0.005099823128319785355949825238269336481254),
1191            Complex64::new(0.4244534840871830045021143490355372016428, 0.002820278933186814021399602648373095266538),
1192            Complex64::new(-0.1021340733271046543881236523269967674156, -0.00001045696456072005761498961861088944159916),
1193            Complex64::new(-0.01000200120119206748855061636187197886859, 0.9805885888237419500266621041508714123763e-8),
1194            Complex64::new(0.001000002000012000023960527532953151819595, -0.9800058800588007290937355024646722133204e-11),
1195            Complex64::new(0.4244549085628511778373438768121222815752, 0.002935393851311701428647152230552122898291),
1196            Complex64::new(-0.1021340732357117208743299813648493928105, -0.00001088377943049851799938998805451564893540),
1197            Complex64::new(-0.01000200120119126652710792390331206563616, 0.1020612612857282306892368985525393707486e-7),
1198            Complex64::new(0.1000000000007333333333344266666666664457e-5, 0.2000000000001333333333323199999999978819e-5),
1199            Complex64::new(0.1999999999994666666666675199999999990248e-5, 0.0),
1200            Complex64::new(0.3013403889237919660346644392864226952119, 0.0),
1201            Complex64::new(0.02503136792640367194699495234782353186858, 0.0),
1202            Complex64::new(0.002500031251171948248596912483183760683918, 0.0),
1203            Complex64::new(0.0, 0.004900078433419939164774792850907128053308),
1204            Complex64::new(0.0, -0.005100088434920074173454208832365950009419),
1205            Complex64::new(0.0, 0.2000000000005333333333341866666666676419e-5),
1206            Complex64::new(0.0, -48.16001211429122974789822893525016528191),
1207            Complex64::new(0.0, 0.4627407029504443513654142715903005954668e174),
1208            Complex64::new(0.0, f64::NEG_INFINITY),
1209            Complex64::new(0.0, 0.0),
1210            Complex64::new(-0.0, 0.0),
1211            Complex64::new(0.0, f64::INFINITY),
1212            Complex64::new(0.0, f64::NEG_INFINITY),
1213            Complex64::new(f64::NAN, f64::NAN),
1214            Complex64::new(f64::NAN, f64::NAN),
1215            Complex64::new(f64::NAN, f64::NAN),
1216            Complex64::new(f64::NAN, 0.0),
1217            Complex64::new(0.0, f64::NAN),
1218            Complex64::new(f64::NAN, f64::NAN),
1219            Complex64::new(f64::NAN, f64::NAN),
1220            Complex64::new(0.01282473148489433743567240624939698290584, -0.2105957276516618621447832572909153498104e-7),
1221            Complex64::new(0.01219875253423634378984109995893708152885, -0.1813040560401824664088425926165834355953e-7),
1222            Complex64::new(0.1020408163265306334945473399689037886997e-7, -0.1041232819658476285651490827866174985330e-25),
1223            Complex64::new(0.9803921568627452865036825956835185367356e-8, -0.9227220299884665067601095648451913375754e-26),
1224            Complex64::new(0.5000000000000000002500000000000000003750e-9, -0.1200000000000000001800000188712838420241e-29),
1225            Complex64::new(5.00000000000000000000025000000000000000000003e-12, -1.20000000000000000000018000000000000000000004e-36),
1226            Complex64::new(5.00000000000000000000000002500000000000000000e-14, -1.20000000000000000000000001800000000000000000e-42),
1227            Complex64::new(5.0e-301, 0.0) 
1228        ];
1229
1230        let tolerance = 1.0e-13;
1231        for n in 0..z.len() {
1232            let computed_dawson = dawson_with_relerror(z[n], 0.);
1233            let relative_error_re = relerr(computed_dawson.re, expected_dawson[n].re);
1234            let relative_error_im = relerr(computed_dawson.im, expected_dawson[n].im);
1235            assert!(relative_error_re < tolerance);
1236            assert!(relative_error_im < tolerance);
1237        }
1238    }
1239
1240} 
1241
errorfunctions/complexerrorfunctions.rs

errorfunctions/
complexerrorfunctions.rs
