libcerf.err_fcts source code

1 /**
2 Computate Dawson, Voigt, and several error functions,
3 based on erfcx, im_w_of_x, w_of_z as implemented in separate files.
4 
5 Copyright:
6     © 2012 Massachusetts Institute of Technology,    
7     © 2013 Forschungszentrum Jülich GmbH,
8     © 2014 Ilya Yaroshenko
9  
10 Authors:
11     Steven G. Johnson, core author;
12     Joachim Wuttke, C package maintainer;
13     $(LINK2 http://plus.google.com/+IlyaYaroshenko, Ilya Yaroshenko), D package maintainer
14 
15 License: 
16     Subject to the terms of the MIT license, as written in the included LICENSE.txt file.
17  */
18 module libcerf.err_fcts;
19 
20 import std.complex;
21 import std.math : isNaN, copysign, SQRT2, M_1_PI;
22 import core.stdc.math : sin, cos, erf, fabs, exp, erfc;
23 
24 import libcerf.w_of_z;
25 import libcerf.im_w_of_x;
26 import libcerf.erfcx_;
27 
28 version(unittest)
29 {
30     import libcerf.testutils;
31 }
32 
33 private enum double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
34 private enum double s2pi = 2.5066282746310005024157652848110; // sqrt(2*pi)
35 
36 /**
37 Compute erfcx(z) = exp(z^2) erfc(z),
38 the complex underflow-compensated complementary error function,
39 trivially related to Faddeeva's w_of_z.
40 */
41 Complex!double cerfcx(Complex!double z) @safe @nogc nothrow
42 {
43     return w_of_z(Complex!double(-z.im, z.re));
44 }
45 
46 unittest
47 {
48     static immutable Z = [ Complex!double(1.234,0.5678) ];
49     static immutable W = [ 
50         Complex!double(0.3382187479799972294747793561190487832579,
51           -0.1116077470811648467464927471872945833154)
52     ];
53     commonTest!cerfcx(Z, W);
54 }
55 
56 /**
57 Compute erfi(z) = -i erf(iz), the rotated complex error function.
58 */
59 Complex!double cerfi(Complex!double z) @safe @nogc nothrow
60 {
61     Complex!double e = cerf(Complex!double(-z.im, z.re));
62     return Complex!double(e.im, -e.re);
63 }
64 
65 unittest
66 {
67     static immutable Z = [ Complex!double(1.234,0.5678) ];
68     static immutable W = [ 
69         Complex!double(1.081032284405373149432716643834106923212,
70           1.926775520840916645838949402886591180834)
71     ];
72     commonTest!cerfi(Z, W);
73 }
74 
75 /**
76 Compute erfi(x) = -i erf(ix), the imaginary error function.
77 */
78 double erfi(double x) @safe @nogc nothrow
79 {
80     immutable x2 = x*x;
81     return x2 > 720 ? copysign(double.infinity, x) : exp(x2) * im_w_of_x(x);
82 }
83 
84 unittest {
85     double[2][] tests = [
86         [0.0, 0],
87         [1.0, 1.65042575879754287602533772956136244389],
88         [-1.0, -1.65042575879754287602533772956136244389],
89         [5.0, 8.29827388067680351614622319074691995187273573754139844564e+9],
90         [-5.0, -8.29827388067680351614622319074691995187273573754139844564e+9],
91     ];
92     foreach(test; tests)
93     {
94         auto x = test[0];
95         auto t = erfi(x);
96         auto y = test[1];
97         assert(relativeErrorCheck(t, y));
98     }
99 }
100 
101 /**
102 Compute dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x),
103 Dawson's integral for a real argument.
104 */
105 double dawson(double x) @safe @nogc nothrow
106 {
107     return spi2 * im_w_of_x(x);
108 }
109 
110 /**
111 Compute Voigt's convolution of a Gaussian
112 G(x,sigma) = 1/sqrt(2*pi)/|sigma| * exp(-x^2/2/sigma^2)
113 and a Lorentzian
114 L(x,gamma) = |gamma| / pi / ( x^2 + gamma^2 ),
115 namely
116 voigt(x,sigma,gamma) =
117       \int_{-infty}^{infty} dx' G(x',sigma) L(x-x',gamma)
118 using the relation
119 voigt(x,sigma,gamma) = Re{ w(z) } / sqrt(2*pi) / |sigma|
120 with
121 z = (x+i*|gamma|) / sqrt(2) / |sigma|.
122 
123 Reference: Abramowitz&Stegun (1964), formula (7.4.13).
124 */
125 //Joachim Wuttke, January 2013.
126 double voigt( double x, double sigma, double gamma ) @safe @nogc nothrow
127 {
128 
129     immutable gam = fabs(gamma);
130     immutable sig = fabs(sigma);
131 
132     if (!gam) {
133         if (!sig) {
134             // It's kind of a delta function
135             return x ? 0 : double.infinity;
136         } else {
137             // It's a pure Gaussian
138             return exp( -x*x/2/(sig*sig) ) / s2pi / sig;
139         }
140     } else {
141         if (!sig) {
142             // It's a pure Lorentzian
143             return gam * M_1_PI / (x*x + gam*gam);
144         } else {
145             // Regular case, both parameters are nonzero
146             Complex!double z = Complex!double(x,gam) / double(SQRT2) / sig;
147             return  w_of_z(z).re / s2pi / sig;
148         }
149     }
150 }
151 
152 /**
153 Compute erf(z), the complex error function,
154 using w_of_z except for certain regions.
155 */
156 //Steven G. Johnson, October 2012.
157 Complex!double cerf(Complex!double z) @safe @nogc nothrow
158 {
159     immutable x = z.re, y = z.im;
160 
161     if (!y)
162         return Complex!double(erf(x), y); // preserve sign of 0
163     if (!x) // handle separately for speed & handling of y = Inf or NaN
164         return Complex!double(x, // preserve sign of 0
165                  /* handle y -> Inf limit manually, since
166                     exp(y^2) -> Inf but Im[w(y)] -> 0, so
167                     IEEE will give us a NaN when it should be Inf */
168                  y*y > 720 ? copysign(double.infinity, y)
169                  : exp(y*y) * im_w_of_x(y));
170   
171     immutable mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
172     immutable mIm_z2 = -2*x*y; // Im(-z^2)
173     if (mRe_z2 < -750) // underflow
174         return Complex!double(copysign(1.0, x), 0.0);
175 
176     /* Handle positive and negative x via different formulas,
177        using the mirror symmetries of w, to avoid overflow/underflow
178        problems from multiplying exponentially large and small quantities. */
179     if (x >= 0) {
180         if (x < 8e-2) {
181             if (fabs(y) < 1e-2)
182                 goto taylor;
183             else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
184                 goto taylor_erfi;
185         }
186         /* don't use complex exp function, since that will produce spurious NaN
187            values when multiplying w in an overflow situation. */
188         return 1.0 - exp(mRe_z2) *
189             (Complex!double(cos(mIm_z2), sin(mIm_z2))
190              * w_of_z(Complex!double(-y,x)));
191     }
192     else { // x < 0
193         if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
194             if (fabs(y) < 1e-2)
195                 goto taylor;
196             else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
197                 goto taylor_erfi;
198         }
199         else if (isNaN(x))
200             return Complex!double(double.nan, y == 0 ? 0 : double.nan);
201         /* don't use complex exp function, since that will produce spurious NaN
202            values when multiplying w in an overflow situation. */
203         return exp(mRe_z2) *
204             (Complex!double(cos(mIm_z2), sin(mIm_z2))
205              * w_of_z(Complex!double(y,-x))) - 1.0;
206     }
207 
208     // Use Taylor series for small |z|, to avoid cancellation inaccuracy
209     //   erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
210 taylor:
211     Complex!double mz2 = Complex!double(mRe_z2, mIm_z2); // -z^2
212     return z * (1.1283791670955125739
213                 + mz2 * (0.37612638903183752464
214                          + mz2 * (0.11283791670955125739
215                                   + mz2 * (0.026866170645131251760
216                                            + mz2 * 0.0052239776254421878422))));
217 
218     /* for small |x| and small |xy|, 
219        use Taylor series to avoid cancellation inaccuracy:
220        erf(x+iy) = erf(iy)
221        + 2*exp(y^2)/sqrt(pi) *
222        [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 
223        - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
224        where:
225        erf(iy) = exp(y^2) * Im[w(y)]
226     */
227 taylor_erfi:
228     double x2 = x*x, y2 = y*y;
229     double expy2 = exp(y2);
230     return Complex!double
231         (expy2 * x * (1.1283791670955125739
232                       - x2 * (0.37612638903183752464
233                               + 0.75225277806367504925*y2)
234                       + x2*x2 * (0.11283791670955125739
235                                  + y2 * (0.45135166683820502956
236                                          + 0.15045055561273500986*y2))),
237          expy2 * (im_w_of_x(y)
238                   - x2*y * (1.1283791670955125739 
239                             - x2 * (0.56418958354775628695 
240                                     + 0.37612638903183752464*y2))));
241 }
242 
243 unittest
244 {
245     static immutable Z = [
246         Complex!double(1,2),
247         Complex!double(-1,2),
248         Complex!double(1,-2),
249         Complex!double(-1,-2),
250         Complex!double(9,-28),
251         Complex!double(21,-33),
252         Complex!double(1e3,1e3),
253         Complex!double(-3001,-1000),
254         Complex!double(1e160,-1e159),
255         Complex!double(5.1e-3, 1e-8),
256         Complex!double(-4.9e-3, 4.95e-3),
257         Complex!double(4.9e-3, 0.5),
258         Complex!double(4.9e-4, -0.5e1),
259         Complex!double(-4.9e-5, -0.5e2),
260         Complex!double(5.1e-3, 0.5),
261         Complex!double(5.1e-4, -0.5e1),
262         Complex!double(-5.1e-5, -0.5e2),
263         Complex!double(1e-6,2e-6),
264         Complex!double(0,2e-6),
265         Complex!double(0,2),
266         Complex!double(0,20),
267         Complex!double(0,200),
268         Complex!double(double.infinity,0),
269         Complex!double(-double.infinity,0),
270         Complex!double(0,double.infinity),
271         Complex!double(0,-double.infinity),
272         Complex!double(double.infinity,double.infinity),
273         Complex!double(double.infinity,-double.infinity),
274         Complex!double(double.nan,double.nan),
275         Complex!double(double.nan,0),
276         Complex!double(0,double.nan),
277         Complex!double(double.nan,double.infinity),
278         Complex!double(double.infinity,double.nan),
279         Complex!double(1e-3,double.nan),
280         Complex!double(7e-2,7e-2),
281         Complex!double(7e-2,-7e-4),
282         Complex!double(-9e-2,7e-4),
283         Complex!double(-9e-2,9e-2),
284         Complex!double(-7e-4,9e-2),
285         Complex!double(7e-2,0.9e-2),
286         Complex!double(7e-2,1.1e-2)
287     ];
288     static immutable W = [ 
289         Complex!double(-0.5366435657785650339917955593141927494421,
290           -5.049143703447034669543036958614140565553),
291         Complex!double(0.5366435657785650339917955593141927494421,
292           -5.049143703447034669543036958614140565553),
293         Complex!double(-0.5366435657785650339917955593141927494421,
294           5.049143703447034669543036958614140565553),
295         Complex!double(0.5366435657785650339917955593141927494421,
296           5.049143703447034669543036958614140565553),
297         Complex!double(0.3359473673830576996788000505817956637777e304,
298           -0.1999896139679880888755589794455069208455e304),
299         Complex!double(0.3584459971462946066523939204836760283645e278,
300           0.3818954885257184373734213077678011282505e280),
301         Complex!double(0.9996020422657148639102150147542224526887,
302           0.00002801044116908227889681753993542916894856),
303         Complex!double(-1, 0),
304         Complex!double(1, 0),
305         Complex!double(0.005754683859034800134412990541076554934877,
306           0.1128349818335058741511924929801267822634e-7),
307         Complex!double(-0.005529149142341821193633460286828381876955,
308           0.005585388387864706679609092447916333443570),
309         Complex!double(0.007099365669981359632319829148438283865814,
310           0.6149347012854211635026981277569074001219),
311         Complex!double(0.3981176338702323417718189922039863062440e8,
312           -0.8298176341665249121085423917575122140650e10),
313         Complex!double(-double.infinity,
314           -double.infinity),
315         Complex!double(0.007389128308257135427153919483147229573895,
316           0.6149332524601658796226417164791221815139),
317         Complex!double(0.4143671923267934479245651547534414976991e8,
318           -0.8298168216818314211557046346850921446950e10),
319         Complex!double(-double.infinity,
320           -double.infinity),
321         Complex!double(0.1128379167099649964175513742247082845155e-5,
322           0.2256758334191777400570377193451519478895e-5),
323         Complex!double(0,
324           0.2256758334194034158904576117253481476197e-5),
325         Complex!double(0,
326           18.56480241457555259870429191324101719886),
327         Complex!double(0,
328           0.1474797539628786202447733153131835124599e173),
329         Complex!double(0,
330           double.infinity),
331         Complex!double(1,0),
332         Complex!double(-1,0),
333         Complex!double(0,double.infinity),
334         Complex!double(0,-double.infinity),
335         Complex!double(double.nan,double.nan),
336         Complex!double(double.nan,double.nan),
337         Complex!double(double.nan,double.nan),
338         Complex!double(double.nan,0),
339         Complex!double(0,double.nan),
340         Complex!double(double.nan,double.nan),
341         Complex!double(double.nan,double.nan),
342         Complex!double(double.nan,double.nan),
343         Complex!double(0.07924380404615782687930591956705225541145,
344           0.07872776218046681145537914954027729115247),
345         Complex!double(0.07885775828512276968931773651224684454495,
346           -0.0007860046704118224342390725280161272277506),
347         Complex!double(-0.1012806432747198859687963080684978759881,
348           0.0007834934747022035607566216654982820299469),
349         Complex!double(-0.1020998418798097910247132140051062512527,
350           0.1010030778892310851309082083238896270340),
351         Complex!double(-0.0007962891763147907785684591823889484764272,
352           0.1018289385936278171741809237435404896152),
353         Complex!double(0.07886408666470478681566329888615410479530,
354           0.01010604288780868961492224347707949372245),
355         Complex!double(0.07886723099940260286824654364807981336591,
356           0.01235199327873258197931147306290916629654)
357     ];
358     commonTest!cerf(Z, W);
359 }
360 
361 
362 /**
363 Compute erfc(z) = 1 - erf(z), the complex complementary error function,
364 using w_of_z except for certain regions.
365 */
366 //Steven G. Johnson, October 2012.
367 Complex!double cerfc(Complex!double z) @safe @nogc nothrow
368 {
369 
370     immutable x = z.re, y = z.im;
371 
372     if (!x)
373         return Complex!double(1,
374                  /* handle y -> Inf limit manually, since
375                     exp(y^2) -> Inf but Im[w(y)] -> 0, so
376                     IEEE will give us a double.nan when it should be Inf */
377                  y*y > 720 ? -copysign(double.infinity, y)
378                  : -exp(y*y) * im_w_of_x(y));
379     if (!y) {
380         if (x*x > 750) // underflow
381             return Complex!double(x >= 0 ? 0.0 : 2.0,
382                      -y); // preserve sign of 0
383         return Complex!double(x >= 0 ? exp(-x*x) * erfcx(x) 
384                  : 2. - exp(-x*x) * erfcx(-x),
385                  -y); // preserve sign of zero
386     }
387 
388     immutable mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
389     immutable mIm_z2 = -2*x*y; // Im(-z^2)
390     if (mRe_z2 < -750) // underflow
391         return Complex!double(x >= 0 ? 0.0 : 2.0, 0);
392 
393     if (x >= 0)
394         return cexp(Complex!double(mRe_z2, mIm_z2))
395             * w_of_z(Complex!double(-y,x));
396     else
397         return 2.0 - cexp(Complex!double(mRe_z2, mIm_z2))
398             * w_of_z(Complex!double(y,-x));
399 }
400 
401 unittest
402 {
403     static immutable Z = [
404         Complex!double(1,2),
405         Complex!double(-1,2),
406         Complex!double(1,-2),
407         Complex!double(-1,-2),
408         Complex!double(9,-28),
409         Complex!double(21,-33),
410         Complex!double(1e3,1e3),
411         Complex!double(-3001,-1000),
412         Complex!double(1e160,-1e159),
413         Complex!double(5.1e-3, 1e-8),
414         Complex!double(0,2e-6),
415         Complex!double(0,2),
416         Complex!double(0,20),
417         Complex!double(0,200),
418         Complex!double(2e-6,0),
419         Complex!double(2,0),
420         Complex!double(20,0),
421         Complex!double(200,0),
422         Complex!double(double.infinity,0),
423         Complex!double(-double.infinity,0),
424         Complex!double(0,double.infinity),
425         Complex!double(0,-double.infinity),
426         Complex!double(double.infinity,double.infinity),
427         Complex!double(double.infinity,-double.infinity),
428         Complex!double(double.nan,double.nan),
429         Complex!double(double.nan,0),
430         Complex!double(0,double.nan),
431         Complex!double(double.nan,double.infinity),
432         Complex!double(double.infinity,double.nan),
433         Complex!double(88,0)
434     ];
435     static immutable W = [
436         Complex!double(1.536643565778565033991795559314192749442,
437           5.049143703447034669543036958614140565553),
438         Complex!double(0.4633564342214349660082044406858072505579,
439           5.049143703447034669543036958614140565553),
440         Complex!double(1.536643565778565033991795559314192749442,
441           -5.049143703447034669543036958614140565553),
442         Complex!double(0.4633564342214349660082044406858072505579,
443           -5.049143703447034669543036958614140565553),
444         Complex!double(-0.3359473673830576996788000505817956637777e304,
445           0.1999896139679880888755589794455069208455e304),
446         Complex!double(-0.3584459971462946066523939204836760283645e278,
447           -0.3818954885257184373734213077678011282505e280),
448         Complex!double(0.0003979577342851360897849852457775473112748,
449           -0.00002801044116908227889681753993542916894856),
450         Complex!double(2, 0),
451         Complex!double(0, 0),
452         Complex!double(0.9942453161409651998655870094589234450651,
453           -0.1128349818335058741511924929801267822634e-7),
454         Complex!double(1,
455           -0.2256758334194034158904576117253481476197e-5),
456         Complex!double(1,
457           -18.56480241457555259870429191324101719886),
458         Complex!double(1,
459           -0.1474797539628786202447733153131835124599e173),
460         Complex!double(1, -double.infinity),
461         Complex!double(0.9999977432416658119838633199332831406314,
462           0),
463         Complex!double(0.004677734981047265837930743632747071389108,
464           0),
465         Complex!double(0.5395865611607900928934999167905345604088e-175,
466           0),
467         Complex!double(0, 0),
468         Complex!double(0, 0),
469         Complex!double(2, 0),
470         Complex!double(1, -double.infinity),
471         Complex!double(1, double.infinity),
472         Complex!double(double.nan, double.nan),
473         Complex!double(double.nan, double.nan),
474         Complex!double(double.nan, double.nan),
475         Complex!double(double.nan, 0),
476         Complex!double(1, double.nan),
477         Complex!double(double.nan, double.nan),
478         Complex!double(double.nan, double.nan),
479         Complex!double(0,0)
480     ];
481     commonTest!cerfc(Z, W);
482 }
483 
484 /**
485 Compute Dawson(z) = sqrt(pi)/2  *  exp(-z^2) * erfi(z),
486 Dawson's integral for a complex argument,
487 using w_of_z except for certain regions.
488 */
489 //Steven G. Johnson, October 2012.
490 Complex!double cdawson(Complex!double z) @safe @nogc nothrow
491 {
492 
493     immutable x = z.re, y = z.im;
494 
495     // handle axes separately for speed & proper handling of x or y = Inf or double.nan
496     if (!y)
497         return Complex!double(spi2 * im_w_of_x(x), -y); // preserve sign of 0
498     if (!x) {
499         immutable y2 = y*y;
500         if (y2 < 2.5e-5) { // Taylor expansion
501             return Complex!double(x, // preserve sign of 0
502                      y * (1.
503                           + y2 * (0.6666666666666666666666666666666666666667
504                                   + y2 * 0.26666666666666666666666666666666666667)));
505         }
506         return Complex!double(x, // preserve sign of 0
507                  spi2 * (y >= 0 
508                          ? exp(y2) - erfcx(y)
509                          : erfcx(-y) - exp(y2)));
510     }
511 
512     immutable mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
513     immutable mIm_z2 = -2*x*y; // Im(-z^2)
514     Complex!double mz2 = Complex!double(mRe_z2, mIm_z2); // -z^2
515 
516     /* Handle positive and negative x via different formulas,
517        using the mirror symmetries of w, to avoid overflow/underflow
518        problems from multiplying exponentially large and small quantities. */
519     if (y >= 0) {
520         if (fabs(y) < 5e-3) {
521             if (fabs(x) < 5e-3)
522                 goto taylor;
523             else if (fabs(mIm_z2) < 5e-3)
524                 goto taylor_realaxis;
525         }
526         Complex!double res = cexp(mz2) - w_of_z(z);
527         return spi2 * Complex!double(-res.im, res.re);
528     }
529     else { // y < 0
530         if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
531             if (fabs(x) < 5e-3)
532                 goto taylor;
533             else if (fabs(mIm_z2) < 5e-3)
534                 goto taylor_realaxis;
535         }
536         else if (isNaN(y))
537             return Complex!double(!x ? 0 : double.nan, double.nan);
538         Complex!double res = w_of_z(-z) - cexp(mz2);
539         return spi2 * Complex!double(-res.im, res.re);
540     }
541 
542     // Use Taylor series for small |z|, to avoid cancellation inaccuracy
543     //     dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
544 taylor:
545     return z * (1.
546                 + mz2 * (0.6666666666666666666666666666666666666667
547                          + mz2 * 0.2666666666666666666666666666666666666667));
548 
549     /* for small |y| and small |xy|, 
550        use Taylor series to avoid cancellation inaccuracy:
551        dawson(x + iy)
552        = D + y^2 (D + x - 2Dx^2)
553        + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
554        + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
555        + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
556        where D = dawson(x) 
557 
558        However, for large |x|, 2Dx -> 1 which gives cancellation problems in
559        this series (many of the leading terms cancel).  So, for large |x|,
560        we need to substitute a continued-fraction expansion for D.
561 
562        dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))
563 
564        The 6 terms shown here seems to be the minimum needed to be
565        accurate as soon as the simpler Taylor expansion above starts
566        breaking down.  Using this 6-term expansion, factoring out the
567        denominator, and simplifying with Maple, we obtain:
568 
569        Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x
570        = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
571        Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y
572        = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4
573 
574        Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction
575        expansion for the real part, and a 2-term expansion for the imaginary
576        part.  (This avoids overflow problems for huge |x|.)  This yields:
577      
578        Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
579        Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)
580 
581     */
582 taylor_realaxis:
583     immutable x2 = x*x;
584     if (x2 > 1600) // |x| > 40
585     { 
586         immutable y2 = y*y;
587         if (x2 > 25e14) {// |x| > 5e7
588             immutable xy2 = (x*y)*(x*y);
589             return Complex!double((0.5 + y2 * (0.5 + 0.25*y2
590                                   - 0.16666666666666666667*xy2)) / x,
591                      y * (-1 + y2 * (-0.66666666666666666667
592                                      + 0.13333333333333333333*xy2
593                                      - 0.26666666666666666667*y2))
594                      / (2*x2 - 1));
595         }
596         return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
597             Complex!double(x * (33 + x2 * (-28 + 4*x2)
598                    + y2 * (18 - 4*x2 + 4*y2)),
599               y * (-15 + x2 * (24 - 4*x2)
600                    + y2 * (4*x2 - 10 - 4*y2)));
601     }
602     else 
603     {
604         immutable D = spi2 * im_w_of_x(x);
605         immutable y2 = y*y;
606         return Complex!double
607             (D + y2 * (D + x - 2*D*x2)
608              + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
609                         + x * (0.83333333333333333333
610                                - 0.33333333333333333333 * x2)),
611              y * (1 - 2*D*x
612                   + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
613                   + y2*y2 * (0.26666666666666666667 -
614                              x2 * (0.6 - 0.13333333333333333333 * x2)
615                              - D*x * (1 - x2 * (1.3333333333333333333
616                                                 - 0.26666666666666666667 * x2)))));
617     }
618 }
619 
620 unittest
621 {
622     static immutable Z = [
623         Complex!double(2,1),
624         Complex!double(-2,1),
625         Complex!double(2,-1),
626         Complex!double(-2,-1),
627         Complex!double(-28,9),
628         Complex!double(33,-21),
629         Complex!double(1e3,1e3),
630         Complex!double(-1000,-3001),
631         Complex!double(1e-8, 5.1e-3),
632         Complex!double(4.95e-3, -4.9e-3),
633         Complex!double(5.1e-3, 5.1e-3),
634         Complex!double(0.5, 4.9e-3),
635         Complex!double(-0.5e1, 4.9e-4),
636         Complex!double(-0.5e2, -4.9e-5),
637         Complex!double(0.5e3, 4.9e-6),
638         Complex!double(0.5, 5.1e-3),
639         Complex!double(-0.5e1, 5.1e-4),
640         Complex!double(-0.5e2, -5.1e-5),
641         Complex!double(1e-6,2e-6),
642         Complex!double(2e-6,0),
643         Complex!double(2,0),
644         Complex!double(20,0),
645         Complex!double(200,0),
646         Complex!double(0,4.9e-3),
647         Complex!double(0,-5.1e-3),
648         Complex!double(0,2e-6),
649         Complex!double(0,-2),
650         Complex!double(0,20),
651         Complex!double(0,-200),
652         Complex!double(double.infinity,0),
653         Complex!double(-double.infinity,0),
654         Complex!double(0,double.infinity),
655         Complex!double(0,-double.infinity),
656         Complex!double(double.infinity,double.infinity),
657         Complex!double(double.infinity,-double.infinity),
658         Complex!double(double.nan,double.nan),
659         Complex!double(double.nan,0),
660         Complex!double(0,double.nan),
661         Complex!double(double.nan,double.infinity),
662         Complex!double(double.infinity,double.nan),
663         Complex!double(39, 6.4e-5),
664         Complex!double(41, 6.09e-5),
665         Complex!double(4.9e7, 5e-11),
666         Complex!double(5.1e7, 4.8e-11),
667         Complex!double(1e9, 2.4e-12),
668         Complex!double(1e11, 2.4e-14),
669         Complex!double(1e13, 2.4e-16),
670         Complex!double(1e300, 2.4e-303)
671     ];
672     static immutable W = [
673         Complex!double(0.1635394094345355614904345232875688576839,
674           -0.1531245755371229803585918112683241066853),
675         Complex!double(-0.1635394094345355614904345232875688576839,
676           -0.1531245755371229803585918112683241066853),
677         Complex!double(0.1635394094345355614904345232875688576839,
678           0.1531245755371229803585918112683241066853),
679         Complex!double(-0.1635394094345355614904345232875688576839,
680           0.1531245755371229803585918112683241066853),
681         Complex!double(-0.01619082256681596362895875232699626384420,
682           -0.005210224203359059109181555401330902819419),
683         Complex!double(0.01078377080978103125464543240346760257008,
684           0.006866888783433775382193630944275682670599),
685         Complex!double(-0.5808616819196736225612296471081337245459,
686           0.6688593905505562263387760667171706325749),
687         Complex!double(double.infinity,
688           -double.infinity),
689         Complex!double(0.1000052020902036118082966385855563526705e-7,
690           0.005100088434920073153418834680320146441685),
691         Complex!double(0.004950156837581592745389973960217444687524,
692           -0.004899838305155226382584756154100963570500),
693         Complex!double(0.005100176864319675957314822982399286703798,
694           0.005099823128319785355949825238269336481254),
695         Complex!double(0.4244534840871830045021143490355372016428,
696           0.002820278933186814021399602648373095266538),
697         Complex!double(-0.1021340733271046543881236523269967674156,
698           -0.00001045696456072005761498961861088944159916),
699         Complex!double(-0.01000200120119206748855061636187197886859,
700           0.9805885888237419500266621041508714123763e-8),
701         Complex!double(0.001000002000012000023960527532953151819595,
702           -0.9800058800588007290937355024646722133204e-11),
703         Complex!double(0.4244549085628511778373438768121222815752,
704           0.002935393851311701428647152230552122898291),
705         Complex!double(-0.1021340732357117208743299813648493928105,
706           -0.00001088377943049851799938998805451564893540),
707         Complex!double(-0.01000200120119126652710792390331206563616,
708           0.1020612612857282306892368985525393707486e-7),
709         Complex!double(0.1000000000007333333333344266666666664457e-5,
710           0.2000000000001333333333323199999999978819e-5),
711         Complex!double(0.1999999999994666666666675199999999990248e-5,
712           0),
713         Complex!double(0.3013403889237919660346644392864226952119,
714           0),
715         Complex!double(0.02503136792640367194699495234782353186858,
716           0),
717         Complex!double(0.002500031251171948248596912483183760683918,
718           0),
719         Complex!double(0,0.004900078433419939164774792850907128053308),
720         Complex!double(0,-0.005100088434920074173454208832365950009419),
721         Complex!double(0,0.2000000000005333333333341866666666676419e-5),
722         Complex!double(0,-48.16001211429122974789822893525016528191),
723         Complex!double(0,0.4627407029504443513654142715903005954668e174),
724         Complex!double(0,-double.infinity),
725         Complex!double(0,0),
726         Complex!double(-0,0),
727         Complex!double(0, double.infinity),
728         Complex!double(0, -double.infinity),
729         Complex!double(double.nan, double.nan),
730         Complex!double(double.nan, double.nan),
731         Complex!double(double.nan, double.nan),
732         Complex!double(double.nan, 0),
733         Complex!double(0, double.nan),
734         Complex!double(double.nan, double.nan),
735         Complex!double(double.nan, double.nan),
736         Complex!double(0.01282473148489433743567240624939698290584,
737           -0.2105957276516618621447832572909153498104e-7),
738         Complex!double(0.01219875253423634378984109995893708152885,
739           -0.1813040560401824664088425926165834355953e-7),
740         Complex!double(0.1020408163265306334945473399689037886997e-7,
741           -0.1041232819658476285651490827866174985330e-25),
742         Complex!double(0.9803921568627452865036825956835185367356e-8,
743           -0.9227220299884665067601095648451913375754e-26),
744         Complex!double(0.5000000000000000002500000000000000003750e-9,
745           -0.1200000000000000001800000188712838420241e-29),
746         Complex!double(5.00000000000000000000025000000000000000000003e-12,
747           -1.20000000000000000000018000000000000000000004e-36),
748         Complex!double(5.00000000000000000000000002500000000000000000e-14,
749           -1.20000000000000000000000001800000000000000000e-42),
750         Complex!double(5e-301, 0)
751     ];
752     commonTest!cdawson(Z, W);
753 }
754 
755 unittest
756 {
757     specialTest!(cerf, erf, 1e-20);
758     specialTest!(cerfi, erfi, 0);
759     specialTest!(cerfcx, erfcx, 0);
760     specialTest!(cerfc, erfc, 1e-20);
761     specialTest!(cdawson, dawson, 1e-20);
762 }