Mercurial > octave-nkf

changeset 15801:2e30a1aadff2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge new upstream Faddeeva release. * Faddeeva.cc: Use numeric_limits for Inf and NaN, support C11 isnan/isinf, use gnulib floor and copysign, convert tabs to spaces, slight accuracy improvements to erf and dawson functions for |z| around 1e-2.
author Steven G. Johnson <stevenj@alum.mit.edu>
date Sat, 15 Dec 2012 14:12:47 -0500
parents 06832c90ae7d
children b6b95d041813
files liboctave/cruft/Faddeeva/Faddeeva.cc
diffstat 1 files changed, 860 insertions(+), 774 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/cruft/Faddeeva/Faddeeva.cc	Fri Dec 14 16:03:50 2012 -0800
+++ b/liboctave/cruft/Faddeeva/Faddeeva.cc	Sat Dec 15 14:12:47 2012 -0500
@@ -1,3 +1,5 @@
+//  -*- mode:c++; tab-width:2; indent-tabs-mode:nil;  -*-
+
 /* Copyright (c) 2012 Massachusetts Institute of Technology
  * 
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -98,40 +100,45 @@
        4 October 2012: Initial public release (SGJ)
        5 October 2012: Revised (SGJ) to fix spelling error,
                        start summation for large x at round(x/a) (> 1)
-		       rather than ceil(x/a) as in the original
-		       paper, which should slightly improve performance
-     		       (and, apparently, slightly improves accuracy)
+                       rather than ceil(x/a) as in the original
+                       paper, which should slightly improve performance
+                       (and, apparently, slightly improves accuracy)
       19 October 2012: Revised (SGJ) to fix bugs for large x, large -y,
                        and 15<x<26. Performance improvements. Prototype
-		       now supplies default value for relerr.
+                       now supplies default value for relerr.
       24 October 2012: Switch to continued-fraction expansion for
                        sufficiently large z, for performance reasons.
-		       Also, avoid spurious overflow for |z| > 1e154.
-		       Set relerr argument to min(relerr,0.1).
+                       Also, avoid spurious overflow for |z| > 1e154.
+                       Set relerr argument to min(relerr,0.1).
       27 October 2012: Enhance accuracy in Re[w(z)] taken by itself,
                        by switching to Alg. 916 in a region near
-		       the real-z axis where continued fractions
-		       have poor relative accuracy in Re[w(z)].  Thanks
-		       to M. Zaghloul for the tip.
+                       the real-z axis where continued fractions
+                       have poor relative accuracy in Re[w(z)].  Thanks
+                       to M. Zaghloul for the tip.
       29 October 2012: Replace SLATEC-derived erfcx routine with
                        completely rewritten code by me, using a very
-		       different algorithm which is much faster.
+                       different algorithm which is much faster.
       30 October 2012: Implemented special-case code for real z
                        (where real part is exp(-x^2) and imag part is
-		        Dawson integral), using algorithm similar to erfx.
-		       Export ImFaddeeva_w function to make Dawson's
-		       integral directly accessible.
+                        Dawson integral), using algorithm similar to erfx.
+                       Export ImFaddeeva_w function to make Dawson's
+                       integral directly accessible.
       3 November 2012: Provide implementations of erf, erfc, erfcx,
                        and Dawson functions in Faddeeva:: namespace,
-		       in addition to Faddeeva::w.  Provide header
-		       file Faddeeva.hh.
+                       in addition to Faddeeva::w.  Provide header
+                       file Faddeeva.hh.
       4 November 2012: Slightly faster erf for real arguments.
                        Updated MATLAB and Octave plugins.
      27 November 2012: Support compilation with either C++ or
                        plain C (using C99 complex numbers).
-		       For real x, use standard-library erf(x)
-		       and erfc(x) if available (for C99 or C++11).
-		       #include "config.h" if HAVE_CONFIG_H is #defined.
+                       For real x, use standard-library erf(x)
+                       and erfc(x) if available (for C99 or C++11).
+                       #include "config.h" if HAVE_CONFIG_H is #defined.
+     15 December 2012: Portability fixes (copysign, Inf/NaN creation),
+                       use CMPLX/__builtin_complex if available in C,
+                       slight accuracy improvements to erf and dawson
+                       functions near the origin.  Use gnulib functions
+                       if GNULIB_NAMESPACE is defined.
 */
 
 /////////////////////////////////////////////////////////////////////////
@@ -145,19 +152,21 @@
 #endif
 
 /////////////////////////////////////////////////////////////////////////
-// basic math stuff, C++ & C99 compatibility
-
-#define Inf (1./0.) // infinity
-#define NaN (0./0.) // NaN
+// macros to allow us to use either C++ or C (with C99 features)
 
 #ifdef __cplusplus
 
 #  include "Faddeeva.hh"
 
-# include <cfloat>
-# include <cmath>
+#  include <cfloat>
+#  include <cmath>
+#  include <limits>
 using namespace std;
 
+// use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS)
+#  define Inf numeric_limits<double>::infinity()
+#  define NaN numeric_limits<double>::quiet_NaN()
+
 typedef complex<double> cmplx;
 
 // Use C-like complex syntax, since the C syntax is more restrictive
@@ -166,7 +175,7 @@
 #  define cimag(z) imag(z)
 #  define cpolar(r,t) polar(r,t)
 
-#  define CMPLX(a,b) cmplx(a,b)
+#  define C(a,b) cmplx(a,b)
 
 #  define FADDEEVA(name) Faddeeva::name
 #  define FADDEEVA_RE(name) Faddeeva::name
@@ -177,31 +186,84 @@
 #    define isnan my_isnan
 static inline bool my_isinf(double x) { return 1/x == 0.; }
 #    define isinf my_isinf
+#  elif (__cplusplus >= 201103L)
+// g++ gets confused between the C and C++ isnan/isinf functions
+#    define isnan std::isnan
+#    define isinf std::isinf
+#  endif
+
+// copysign was introduced in C++11 (and is also in POSIX and C99)
+#  if defined(_WIN32) || defined(__WIN32__)
+#    define copysign _copysign // of course MS had to be different
+#  elif defined(GNULIB_NAMESPACE) // we are using using gnulib <cmath>
+#    define copysign GNULIB_NAMESPACE::copysign
+#  elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX)
+static inline double my_copysign(double x, double y) { return y<0 ? -x : x; }
+#    define copysign my_copysign
+#  endif
+
+// If we are using the gnulib <cmath> (e.g. in the GNU Octave sources),
+// gnulib generates a link warning if we use ::floor instead of gnulib::floor.
+// This warning is completely innocuous because the only difference between
+// gnulib::floor and the system ::floor (and only on ancient OSF systems)
+// has to do with floor(-0), which doesn't occur in the usage below, but
+// the Octave developers prefer that we silence the warning.
+#  ifdef GNULIB_NAMESPACE
+#    define floor GNULIB_NAMESPACE::floor
 #  endif
 
 #else // !__cplusplus, i.e. pure C (requires C99 features)
 
 #  include "Faddeeva.h"
 
+#  define _GNU_SOURCE // enable GNU libc NAN extension if possible
+
 #  include <float.h>
 #  include <math.h>
 
 typedef double complex cmplx;
 
-#  define CMPLX(a,b) ((a) + I*(b))
+#  define FADDEEVA(name) Faddeeva_ ## name
+#  define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
+
+/* Constructing complex numbers like 0+i*NaN is problematic in C99
+   without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
+   I is a complex (rather than imaginary) constant.  For some reason,
+   however, it works fine in (pre-4.7) gcc if I define Inf and NaN as
+   1/0 and 0/0 (and only if I compile with optimization -O1 or more),
+   but not if I use the INFINITY or NAN macros. */
+
+/* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro
+   may not be defined unless we are using a recent (2012) version of
+   glibc and compile with -std=c11... note that icc lies about being
+   gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
+#  if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER))
+#    define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
+#  endif
+
+#  ifdef CMPLX // C11
+#    define C(a,b) CMPLX(a,b)
+#    define Inf INFINITY // C99 infinity
+#    ifdef NAN // GNU libc extension
+#      define NaN NAN
+#    else
+#      define NaN (0./0.) // NaN
+#    endif
+#  else
+#    define C(a,b) ((a) + I*(b))
+#    define Inf (1./0.) 
+#    define NaN (0./0.) 
+#  endif
 
 static inline cmplx cpolar(double r, double t)
 {
   if (r == 0.0 && !isnan(t))
     return 0.0;
   else
-    return CMPLX(r * cos(t), r * sin(t));
+    return C(r * cos(t), r * sin(t));
 }
 
-#  define FADDEEVA(name) Faddeeva_ ## name
-#  define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
-
-#endif
+#endif // !__cplusplus, i.e. pure C (requires C99 features)
 
 /////////////////////////////////////////////////////////////////////////
 // Auxiliary routines to compute other special functions based on w(z)
@@ -209,7 +271,7 @@
 // compute erfcx(z) = exp(z^2) erfz(z)
 cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
 {
-  return FADDEEVA(w)(CMPLX(-cimag(z), creal(z)), relerr);
+  return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
 }
 
 // compute the error function erf(x)
@@ -225,20 +287,22 @@
     return (x >= 0 ? 1.0 : -1.0);
 
   if (x >= 0) {
-    if (x < 5e-3) goto taylor;
+    if (x < 8e-2) goto taylor;
     return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
   }
   else { // x < 0
-    if (x > -5e-3) goto taylor;
+    if (x > -8e-2) goto taylor;
     return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0;
   }
 
   // Use Taylor series for small |x|, to avoid cancellation inaccuracy
-  //     erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - ...)
+  //   erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...)
  taylor:
   return x * (1.1283791670955125739
-	      + mx2 * (0.37612638903183752464
-		       + mx2 * 0.11283791670955125739));
+              + mx2 * (0.37612638903183752464
+                       + mx2 * (0.11283791670955125739
+                                + mx2 * (0.026866170645131251760
+                                         + mx2 * 0.0052239776254421878422))));
 #endif
 }
 
@@ -248,16 +312,16 @@
   double x = creal(z), y = cimag(z);
 
   if (y == 0)
-    return CMPLX(FADDEEVA_RE(erf)(x),
-			   y); // preserve sign of 0
+    return C(FADDEEVA_RE(erf)(x),
+             y); // preserve sign of 0
   if (x == 0) // handle separately for speed & handling of y = Inf or NaN
-    return CMPLX(x, // preserve sign of 0
-			   /* handle y -> Inf limit manually, since
-			      exp(y^2) -> Inf but Im[w(y)] -> 0, so
-			      IEEE will give us a NaN when it should be Inf */
-			   y*y > 720 ? (y > 0 ? Inf : -Inf)
-			   : exp(y*y) * FADDEEVA(w_im)(y));
-
+    return C(x, // preserve sign of 0
+             /* handle y -> Inf limit manually, since
+                exp(y^2) -> Inf but Im[w(y)] -> 0, so
+                IEEE will give us a NaN when it should be Inf */
+             y*y > 720 ? (y > 0 ? Inf : -Inf)
+             : exp(y*y) * FADDEEVA(w_im)(y));
+  
   double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
   double mIm_z2 = -2*x*y; // Im(-z^2)
   if (mRe_z2 < -750) // underflow
@@ -267,49 +331,51 @@
      using the mirror symmetries of w, to avoid overflow/underflow
      problems from multiplying exponentially large and small quantities. */
   if (x >= 0) {
-    if (x < 5e-3) {
-      if (fabs(y) < 5e-3)
-	goto taylor;
-      else if (fabs(mIm_z2) < 5e-3)
-	goto taylor_erfi;
+    if (x < 8e-2) {
+      if (fabs(y) < 1e-2)
+        goto taylor;
+      else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
+        goto taylor_erfi;
     }
     /* don't use complex exp function, since that will produce spurious NaN
        values when multiplying w in an overflow situation. */
     return 1.0 - exp(mRe_z2) *
-      (CMPLX(cos(mIm_z2), sin(mIm_z2))
-       * FADDEEVA(w)(CMPLX(-y,x), relerr));
+      (C(cos(mIm_z2), sin(mIm_z2))
+       * FADDEEVA(w)(C(-y,x), relerr));
   }
   else { // x < 0
-    if (x > -5e-3) { // duplicate from above to avoid fabs(x) call
-      if (fabs(y) < 5e-3)
-	goto taylor;
-      else if (fabs(mIm_z2) < 5e-3)
-	goto taylor_erfi;
+    if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
+      if (fabs(y) < 1e-2)
+        goto taylor;
+      else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
+        goto taylor_erfi;
     }
     else if (isnan(x))
-      return CMPLX(NaN, y == 0 ? 0 : NaN);
+      return C(NaN, y == 0 ? 0 : NaN);
     /* don't use complex exp function, since that will produce spurious NaN
        values when multiplying w in an overflow situation. */
     return exp(mRe_z2) *
-      (CMPLX(cos(mIm_z2), sin(mIm_z2))
-       * FADDEEVA(w)(CMPLX(y,-x), relerr)) - 1.0;
+      (C(cos(mIm_z2), sin(mIm_z2))
+       * FADDEEVA(w)(C(y,-x), relerr)) - 1.0;
   }
 
   // Use Taylor series for small |z|, to avoid cancellation inaccuracy
-  //     erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - ...)
+  //   erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
  taylor:
   {
-    cmplx mz2 = CMPLX(mRe_z2, mIm_z2); // -z^2
+    cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
     return z * (1.1283791670955125739
-		+ mz2 * (0.37612638903183752464
-			 + mz2 * 0.11283791670955125739));
+                + mz2 * (0.37612638903183752464
+                         + mz2 * (0.11283791670955125739
+                                  + mz2 * (0.026866170645131251760
+                                          + mz2 * 0.0052239776254421878422))));
   }
 
   /* for small |x| and small |xy|, 
      use Taylor series to avoid cancellation inaccuracy:
        erf(x+iy) = erf(iy)
           + 2*exp(y^2)/sqrt(pi) *
-	    [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 
+            [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 
               - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
      where:
         erf(iy) = exp(y^2) * Im[w(y)]
@@ -318,25 +384,25 @@
   {
     double x2 = x*x, y2 = y*y;
     double expy2 = exp(y2);
-    return CMPLX
+    return C
       (expy2 * x * (1.1283791670955125739
-		    - x2 * (0.37612638903183752464
-			    + 0.75225277806367504925*y2)
-		    + x2*x2 * (0.11283791670955125739
-			       + y2 * (0.45135166683820502956
-				       + 0.15045055561273500986*y2))),
+                    - x2 * (0.37612638903183752464
+                            + 0.75225277806367504925*y2)
+                    + x2*x2 * (0.11283791670955125739
+                               + y2 * (0.45135166683820502956
+                                       + 0.15045055561273500986*y2))),
        expy2 * (FADDEEVA(w_im)(y)
-		- x2*y * (1.1283791670955125739 
-			  - x2 * (0.56418958354775628695 
-				  + 0.37612638903183752464*y2))));
+                - x2*y * (1.1283791670955125739 
+                          - x2 * (0.56418958354775628695 
+                                  + 0.37612638903183752464*y2))));
   }
 }
 
 // erfi(z) = -i erf(iz)
 cmplx FADDEEVA(erfi)(cmplx z, double relerr)
 {
-  cmplx e = FADDEEVA(erf)(CMPLX(-cimag(z),creal(z)), relerr);
-  return CMPLX(cimag(e), -creal(e));
+  cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr);
+  return C(cimag(e), -creal(e));
 }
 
 // erfi(x) = -i erf(ix)
@@ -367,19 +433,19 @@
   double x = creal(z), y = cimag(z);
 
   if (x == 0.)
-    return CMPLX(1,
-			   /* handle y -> Inf limit manually, since
-			      exp(y^2) -> Inf but Im[w(y)] -> 0, so
-			      IEEE will give us a NaN when it should be Inf */
-			   y*y > 720 ? (y > 0 ? -Inf : Inf)
-			   : -exp(y*y) * FADDEEVA(w_im)(y));
+    return C(1,
+             /* handle y -> Inf limit manually, since
+                exp(y^2) -> Inf but Im[w(y)] -> 0, so
+                IEEE will give us a NaN when it should be Inf */
+             y*y > 720 ? (y > 0 ? -Inf : Inf)
+             : -exp(y*y) * FADDEEVA(w_im)(y));
   if (y == 0.) {
     if (x*x > 750) // underflow
-      return CMPLX(x >= 0 ? 0.0 : 2.0,
-		   -y); // preserve sign of 0
-    return CMPLX(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 
-		 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x),
-		 -y); // preserve sign of zero
+      return C(x >= 0 ? 0.0 : 2.0,
+               -y); // preserve sign of 0
+    return C(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 
+             : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x),
+             -y); // preserve sign of zero
   }
 
   double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
@@ -388,11 +454,11 @@
     return (x >= 0 ? 0.0 : 2.0);
 
   if (x >= 0)
-    return cexp(CMPLX(mRe_z2, mIm_z2))
-      * FADDEEVA(w)(CMPLX(-y,x), relerr);
+    return cexp(C(mRe_z2, mIm_z2))
+      * FADDEEVA(w)(C(-y,x), relerr);
   else
-    return 2.0 - cexp(CMPLX(mRe_z2, mIm_z2))
-      * FADDEEVA(w)(CMPLX(y,-x), relerr);
+    return 2.0 - cexp(C(mRe_z2, mIm_z2))
+      * FADDEEVA(w)(C(y,-x), relerr);
 }
 
 // compute Dawson(x) = sqrt(pi)/2  *  exp(-x^2) * erfi(x)
@@ -410,25 +476,25 @@
 
   // handle axes separately for speed & proper handling of x or y = Inf or NaN
   if (y == 0)
-    return CMPLX(spi2 * FADDEEVA(w_im)(x),
-			   -y); // preserve sign of 0
+    return C(spi2 * FADDEEVA(w_im)(x),
+             -y); // preserve sign of 0
   if (x == 0) {
     double y2 = y*y;
     if (y2 < 2.5e-5) { // Taylor expansion
-      return CMPLX(x, // preserve sign of 0
-	 y * (1.
-	      + y2 * (0.6666666666666666666666666666666666666667
-		      + y2 * 0.2666666666666666666666666666666666666667)));
+      return C(x, // preserve sign of 0
+               y * (1.
+                    + y2 * (0.6666666666666666666666666666666666666667
+                            + y2 * 0.26666666666666666666666666666666666667)));
     }
-    return CMPLX(x, // preserve sign of 0
-			   spi2 * (y >= 0 
-				   ? exp(y2) - FADDEEVA_RE(erfcx)(y)
-				   : FADDEEVA_RE(erfcx)(-y) - exp(y2)));
+    return C(x, // preserve sign of 0
+             spi2 * (y >= 0 
+                     ? exp(y2) - FADDEEVA_RE(erfcx)(y)
+                     : FADDEEVA_RE(erfcx)(-y) - exp(y2)));
   }
 
   double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
   double mIm_z2 = -2*x*y; // Im(-z^2)
-  cmplx mz2 = CMPLX(mRe_z2, mIm_z2); // -z^2
+  cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
 
   /* Handle positive and negative x via different formulas,
      using the mirror symmetries of w, to avoid overflow/underflow
@@ -436,32 +502,32 @@
   if (y >= 0) {
     if (y < 5e-3) {
       if (fabs(x) < 5e-3)
-	goto taylor;
+        goto taylor;
       else if (fabs(mIm_z2) < 5e-3)
-	goto taylor_realaxis;
+        goto taylor_realaxis;
     }
     cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr);
-    return spi2 * CMPLX(-cimag(res), creal(res));
+    return spi2 * C(-cimag(res), creal(res));
   }
   else { // y < 0
     if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
       if (fabs(x) < 5e-3)
-	goto taylor;
+        goto taylor;
       else if (fabs(mIm_z2) < 5e-3)
-	goto taylor_realaxis;
+        goto taylor_realaxis;
     }
     else if (isnan(y))
-      return CMPLX(x == 0 ? 0 : NaN, NaN);
+      return C(x == 0 ? 0 : NaN, NaN);
     cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2);
-    return spi2 * CMPLX(-cimag(res), creal(res));
+    return spi2 * C(-cimag(res), creal(res));
   }
 
   // Use Taylor series for small |z|, to avoid cancellation inaccuracy
   //     dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
  taylor:
   return z * (1.
-	      + mz2 * (0.6666666666666666666666666666666666666667
-		       + mz2 * 0.2666666666666666666666666666666666666667));
+              + mz2 * (0.6666666666666666666666666666666666666667
+                       + mz2 * 0.2666666666666666666666666666666666666667));
 
   /* for small |y| and small |xy|, 
      use Taylor series to avoid cancellation inaccuracy:
@@ -502,34 +568,34 @@
     if (x2 > 1600) { // |x| > 40
       double y2 = y*y;
       if (x2 > 25e14) {// |x| > 5e7
-	double xy2 = (x*y)*(x*y);
-	return CMPLX((0.5 + y2 * (0.5 + 0.25*y2
-					    - 0.16666666666666666667*xy2)) / x,
-			       y * (-1 + y2 * (-0.66666666666666666667
-					       + 0.13333333333333333333*xy2
-					       - 0.26666666666666666667*y2))
-			       / (2*x2 - 1));
+        double xy2 = (x*y)*(x*y);
+        return C((0.5 + y2 * (0.5 + 0.25*y2
+                              - 0.16666666666666666667*xy2)) / x,
+                 y * (-1 + y2 * (-0.66666666666666666667
+                                 + 0.13333333333333333333*xy2
+                                 - 0.26666666666666666667*y2))
+                 / (2*x2 - 1));
       }
       return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
-	CMPLX(x * (33 + x2 * (-28 + 4*x2)
-			     + y2 * (18 - 4*x2 + 4*y2)),
-			y * (-15 + x2 * (24 - 4*x2)
-			     + y2 * (4*x2 - 10 - 4*y2)));
+        C(x * (33 + x2 * (-28 + 4*x2)
+               + y2 * (18 - 4*x2 + 4*y2)),
+          y * (-15 + x2 * (24 - 4*x2)
+               + y2 * (4*x2 - 10 - 4*y2)));
     }
     else {
       double D = spi2 * FADDEEVA(w_im)(x);
       double x2 = x*x, y2 = y*y;
-      return CMPLX
-	(D + y2 * (D + x - 2*D*x2)
-	 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
-		    + x * (0.83333333333333333333
-			   - 0.33333333333333333333 * x2)),
-	 y * (1 - 2*D*x
-	      + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
-	      + y2*y2 * (0.26666666666666666667 -
-			 x2 * (0.6 - 0.13333333333333333333 * x2)
-			 - D*x * (1 - x2 * (1.3333333333333333333
-					    - 0.26666666666666666667 * x2)))));
+      return C
+        (D + y2 * (D + x - 2*D*x2)
+         + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
+                    + x * (0.83333333333333333333
+                           - 0.33333333333333333333 * x2)),
+         y * (1 - 2*D*x
+              + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
+              + y2*y2 * (0.26666666666666666667 -
+                         x2 * (0.6 - 0.13333333333333333333 * x2)
+                         - D*x * (1 - x2 * (1.3333333333333333333
+                                            - 0.26666666666666666667 * x2)))));
     }
   }
 }
@@ -545,7 +611,7 @@
 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2
 static inline double sinh_taylor(double x) {
   return x * (1 + (x*x) * (0.1666666666666666666667
-			   + 0.00833333333333333333333 * (x*x)));
+                           + 0.00833333333333333333333 * (x*x)));
 }
 
 static inline double sqr(double x) { return x*x; }
@@ -612,11 +678,11 @@
 cmplx FADDEEVA(w)(cmplx z, double relerr)
 {
   if (creal(z) == 0.0)
-    return CMPLX(FADDEEVA_RE(erfcx)(cimag(z)), 
-			   creal(z)); // give correct sign of 0 in cimag(w)
+    return C(FADDEEVA_RE(erfcx)(cimag(z)), 
+             creal(z)); // give correct sign of 0 in cimag(w)
   else if (cimag(z) == 0)
-    return CMPLX(exp(-sqr(creal(z))),
-			   FADDEEVA(w_im)(creal(z)));
+    return C(exp(-sqr(creal(z))),
+             FADDEEVA(w_im)(creal(z)));
 
   double a, a2, c;
   if (relerr <= DBL_EPSILON) {
@@ -643,11 +709,11 @@
 
 #if USE_CONTINUED_FRACTION
   if (ya > 7 || (x > 6  // continued fraction is faster
-		 /* As pointed out by M. Zaghloul, the continued
-		    fraction seems to give a large relative error in
-		    Re w(z) for |x| ~ 6 and small |y|, so use
-		    algorithm 816 in this region: */
-		 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
+                 /* As pointed out by M. Zaghloul, the continued
+                    fraction seems to give a large relative error in
+                    Re w(z) for |x| ~ 6 and small |y|, so use
+                    algorithm 816 in this region: */
+                 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
     
     /* Poppe & Wijers suggest using a number of terms
            nu = 3 + 1442 / (26*rho + 77)
@@ -663,25 +729,25 @@
     double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
     if (x + ya > 4000) { // nu <= 2
       if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z
-	// scale to avoid overflow
-	if (x > ya) {
-	  double yax = ya / xs; 
-	  double denom = ispi / (xs + yax*ya);
-	  ret = CMPLX(denom*yax, denom);
-	}
-	else if (isinf(ya))
-	  return ((isnan(x) || y < 0) 
-		  ? CMPLX(NaN,NaN) : CMPLX(0,0));
-	else {
-	  double xya = xs / ya;
-	  double denom = ispi / (xya*xs + ya);
-	  ret = CMPLX(denom, denom*xya);
-	}
+        // scale to avoid overflow
+        if (x > ya) {
+          double yax = ya / xs; 
+          double denom = ispi / (xs + yax*ya);
+          ret = C(denom*yax, denom);
+        }
+        else if (isinf(ya))
+          return ((isnan(x) || y < 0) 
+                  ? C(NaN,NaN) : C(0,0));
+        else {
+          double xya = xs / ya;
+          double denom = ispi / (xya*xs + ya);
+          ret = C(denom, denom*xya);
+        }
       }
       else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
-	double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
-	double denom = ispi / (dr*dr + di*di);
-	ret = CMPLX(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
+        double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
+        double denom = ispi / (dr*dr + di*di);
+        ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
       }
     }
     else { // compute nu(z) estimate and do general continued fraction
@@ -689,21 +755,21 @@
       double nu = floor(c0 + c1 / (c2*x + c3*ya + c4));
       double wr = xs, wi = ya;
       for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) {
-	// w <- z - nu/w:
-	double denom = nu / (wr*wr + wi*wi);
-	wr = xs - wr * denom;
-	wi = ya + wi * denom;
+        // w <- z - nu/w:
+        double denom = nu / (wr*wr + wi*wi);
+        wr = xs - wr * denom;
+        wi = ya + wi * denom;
       }
       { // w(z) = i/sqrt(pi) / w:
-	double denom = ispi / (wr*wr + wi*wi);
-	ret = CMPLX(denom*wi, denom*wr);
+        double denom = ispi / (wr*wr + wi*wi);
+        ret = C(denom*wi, denom*wr);
       }
     }
     if (y < 0) {
       // 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) 
-      return 2.0*cexp(CMPLX((ya-xs)*(xs+ya), 2*xs*y)) - ret;
+      return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
     }
     else
       return ret;
@@ -716,18 +782,18 @@
     if (x > ya) {
       double yax = ya / xs; 
       double denom = ispi / (xs + yax*ya);
-      ret = CMPLX(denom*yax, denom);
+      ret = C(denom*yax, denom);
     }
     else {
       double xya = xs / ya;
       double denom = ispi / (xya*xs + ya);
-      ret = CMPLX(denom, denom*xya);
+      ret = C(denom, denom*xya);
     }
     if (y < 0) {
       // 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) 
-      return 2.0*cexp(CMPLX((ya-xs)*(xs+ya), 2*xs*y)) - ret;
+      return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
     }
     else
       return ret;
@@ -751,7 +817,7 @@
     double expx2;
 
     if (isnan(y))
-      return CMPLX(y,y);
+      return C(y,y);
     
     /* Somewhat ugly copy-and-paste duplication here, but I see significant
        speedups from using the special-case code with the precomputed
@@ -759,80 +825,80 @@
 
     if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table
       if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
-	const double x2 = x*x;
-	expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
+        const double x2 = x*x;
+        expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
         // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
-	const double ax2 = 1.036642960860171859744*x; // 2*a*x
-	const double exp2ax =
-	  1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
-	const double expm2ax =
-	  1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
-	for (int n = 1; 1; ++n) {
-	  const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
-	  prod2ax *= exp2ax;
-	  prodm2ax *= expm2ax;
-	  sum1 += coef;
-	  sum2 += coef * prodm2ax;
-	  sum3 += coef * prod2ax;
-	  
-	  // really = sum5 - sum4
-	  sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
-	  
-	  // test convergence via sum3
-	  if (coef * prod2ax < relerr * sum3) break;
-	}
+        const double ax2 = 1.036642960860171859744*x; // 2*a*x
+        const double exp2ax =
+          1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
+        const double expm2ax =
+          1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
+        for (int n = 1; 1; ++n) {
+          const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
+          prod2ax *= exp2ax;
+          prodm2ax *= expm2ax;
+          sum1 += coef;
+          sum2 += coef * prodm2ax;
+          sum3 += coef * prod2ax;
+          
+          // really = sum5 - sum4
+          sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
+          
+          // test convergence via sum3
+          if (coef * prod2ax < relerr * sum3) break;
+        }
       }
       else { // x > 5e-4, compute sum4 and sum5 separately
-	expx2 = exp(-x*x);
-	const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
-	for (int n = 1; 1; ++n) {
-	  const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
-	  prod2ax *= exp2ax;
-	  prodm2ax *= expm2ax;
-	  sum1 += coef;
-	  sum2 += coef * prodm2ax;
-	  sum4 += (coef * prodm2ax) * (a*n);
-	  sum3 += coef * prod2ax;
-	  sum5 += (coef * prod2ax) * (a*n);
-	  // test convergence via sum5, since this sum has the slowest decay
-	  if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
-	}
+        expx2 = exp(-x*x);
+        const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
+        for (int n = 1; 1; ++n) {
+          const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
+          prod2ax *= exp2ax;
+          prodm2ax *= expm2ax;
+          sum1 += coef;
+          sum2 += coef * prodm2ax;
+          sum4 += (coef * prodm2ax) * (a*n);
+          sum3 += coef * prod2ax;
+          sum5 += (coef * prod2ax) * (a*n);
+          // test convergence via sum5, since this sum has the slowest decay
+          if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
+        }
       }
     }
     else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly
       const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
       if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
-	const double x2 = x*x;
-	expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
-	for (int n = 1; 1; ++n) {
-	  const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
-	  prod2ax *= exp2ax;
-	  prodm2ax *= expm2ax;
-	  sum1 += coef;
-	  sum2 += coef * prodm2ax;
-	  sum3 += coef * prod2ax;
-	  
-	  // really = sum5 - sum4
-	  sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
-	  
-	  // test convergence via sum3
-	  if (coef * prod2ax < relerr * sum3) break;
-	}
+        const double x2 = x*x;
+        expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
+        for (int n = 1; 1; ++n) {
+          const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
+          prod2ax *= exp2ax;
+          prodm2ax *= expm2ax;
+          sum1 += coef;
+          sum2 += coef * prodm2ax;
+          sum3 += coef * prod2ax;
+          
+          // really = sum5 - sum4
+          sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
+          
+          // test convergence via sum3
+          if (coef * prod2ax < relerr * sum3) break;
+        }
       }
       else { // x > 5e-4, compute sum4 and sum5 separately
-	expx2 = exp(-x*x);
-	for (int n = 1; 1; ++n) {
-	  const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
-	  prod2ax *= exp2ax;
-	  prodm2ax *= expm2ax;
-	  sum1 += coef;
-	  sum2 += coef * prodm2ax;
-	  sum4 += (coef * prodm2ax) * (a*n);
-	  sum3 += coef * prod2ax;
-	  sum5 += (coef * prod2ax) * (a*n);
-	  // test convergence via sum5, since this sum has the slowest decay
-	  if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
-	}
+        expx2 = exp(-x*x);
+        for (int n = 1; 1; ++n) {
+          const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
+          prod2ax *= exp2ax;
+          prodm2ax *= expm2ax;
+          sum1 += coef;
+          sum2 += coef * prodm2ax;
+          sum4 += (coef * prodm2ax) * (a*n);
+          sum3 += coef * prod2ax;
+          sum5 += (coef * prod2ax) * (a*n);
+          // test convergence via sum5, since this sum has the slowest decay
+          if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
+        }
       }
     }
     const double expx2erfcxy = // avoid spurious overflow for large negative y
@@ -841,7 +907,7 @@
     if (y > 5) { // imaginary terms cancel
       const double sinxy = sin(x*y);
       ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
-	+ (c*x*expx2) * sinxy * sinc(x*y, sinxy);
+        + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
     }
     else {
       double xs = creal(z);
@@ -849,25 +915,25 @@
       const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y);
       const double coef1 = expx2erfcxy - c*y*sum1;
       const double coef2 = c*xs*expx2;
-      ret = CMPLX(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
-			    coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
+      ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
+              coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
     }
   }
   else { // x large: only sum3 & sum5 contribute (see above note)    
     if (isnan(x))
-      return CMPLX(x,x);
+      return C(x,x);
     if (isnan(y))
-      return CMPLX(y,y);
+      return C(y,y);
 
 #if USE_CONTINUED_FRACTION
     ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
 #else
     if (y < 0) {
       /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so
-	 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
-	 if y*y - x*x > -36 or so.  So, compute this term just in case.
-	 We also need the -exp(-x*x) term to compute Re[w] accurately
-	 in the case where y is very small. */
+         erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
+         if y*y - x*x > -36 or so.  So, compute this term just in case.
+         We also need the -exp(-x*x) term to compute Re[w] accurately
+         in the case where y is very small. */
       ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
     }
     else
@@ -899,8 +965,8 @@
     }
   }
  finish:
-  return ret + CMPLX((0.5*c)*y*(sum2+sum3), 
-			       (0.5*c)*copysign(sum5-sum4, creal(z)));
+  return ret + C((0.5*c)*y*(sum2+sum3), 
+                 (0.5*c)*copysign(sum5-sum4, creal(z)));
 }
 
 /////////////////////////////////////////////////////////////////////////
@@ -918,14 +984,14 @@
 
       a) It maps x to y = 4 / (4+x) in [0,1].  This simple transformation,
          inspired by a similar transformation in the octave-forge/specfun
-	 erfcx by Soren Hauberg, results in much faster Chebyshev convergence
-	 than other simple transformations I have examined.
+         erfcx by Soren Hauberg, results in much faster Chebyshev convergence
+         than other simple transformations I have examined.
 
       b) Instead of using a single Chebyshev polynomial for the entire
          [0,1] y interval, we break the interval up into 100 equal
-	 subintervals, with a switch/lookup table, and use much lower
-	 degree Chebyshev polynomials in each subinterval. This greatly
-	 improves performance in my tests.
+         subintervals, with a switch/lookup table, and use much lower
+         degree Chebyshev polynomials in each subinterval. This greatly
+         improves performance in my tests.
 
    For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
    with the usual checks for overflow etcetera.
@@ -1357,16 +1423,16 @@
     if (x > 50) { // continued-fraction expansion is faster
       const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
       if (x > 5e7) // 1-term expansion, important to avoid overflow
-	return ispi / x;
+        return ispi / x;
       /* 5-term expansion (rely on compiler for CSE), simplified from:
-	        ispi / (x+0.5/(x+1/(x+1.5/(x+2/x))))  */
+                ispi / (x+0.5/(x+1/(x+1.5/(x+2/x))))  */
       return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
     }
     return erfcx_y100(400/(4+x));
   }
   else
     return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x) 
-				   : 2*exp(x*x) - erfcx_y100(400/(4-x)));
+                                   : 2*exp(x*x) - erfcx_y100(400/(4-x)));
 }
 
 /////////////////////////////////////////////////////////////////////////
@@ -1777,21 +1843,15 @@
       double t = 2*y100 - 193;
       return 0.40889797115352738582e-1 + (-0.60426484889413678200e-2 + (0.28953496450191694606e-4 + (-0.21982952021823718400e-7 + (-0.11044169117553026211e-8 + (0.13344562332430552171e-10 - 0.10091231402844444444e-12 * t) * t) * t) * t) * t) * t;
     }
-    case 97: {
-      double t = 2*y100 - 195;
-      return 0.28920121009594899986e-1 + (-0.59271325915413781788e-2 + (0.28796136372768177423e-4 + (-0.30300382596279568642e-7 + (-0.97688275022802329749e-9 + (0.12179215701512592356e-10 - 0.93380988481777777779e-13 * t) * t) * t) * t) * t) * t;
-    }
-    case 98: {
-      double t = 2*y100 - 197;
-      return 0.17180782722617876655e-1 + (-0.58123419543161127769e-2 + (0.28591841095380959666e-4 + (-0.37642963496443667043e-7 + (-0.86055809047367300024e-9 + (0.11101709356762665578e-10 - 0.86272947493333333334e-13 * t) * t) * t) * t) * t) * t;
-    }
-  case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.010101...)
-      //  (2/sqrt(pi)) * (x - 2/3 x^3  + 4/15 x^5  - 8/105 x^7) 
+  case 97: case 98:
+  case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
+      //  (2/sqrt(pi)) * (x - 2/3 x^3  + 4/15 x^5  - 8/105 x^7 + 16/945 x^9) 
       double x2 = x*x;
       return x * (1.1283791670955125739
-		  - x2 * (0.75225277806367504925
-			  - x2 * (0.30090111122547001970
-				  - x2 * 0.085971746064420005629)));
+                  - x2 * (0.75225277806367504925
+                          - x2 * (0.30090111122547001970
+                                  - x2 * (0.085971746064420005629
+                                          - x2 * 0.016931216931216931217))));
     }
   }
   /* Since 0 <= y100 < 101, this is only reached if x is NaN,
@@ -1805,9 +1865,9 @@
     if (x > 45) { // continued-fraction expansion is faster
       const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
       if (x > 5e7) // 1-term expansion, important to avoid overflow
-	return ispi / x;
+        return ispi / x;
       /* 5-term expansion (rely on compiler for CSE), simplified from:
-	        ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
+                ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
       return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
     }
     return w_im_y100(100/(1+x), x);
@@ -1816,9 +1876,9 @@
     if (x < -45) { // continued-fraction expansion is faster
       const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
       if (x < -5e7) // 1-term expansion, important to avoid overflow
-	return ispi / x;
+        return ispi / x;
       /* 5-term expansion (rely on compiler for CSE), simplified from:
-	        ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
+                ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
       return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
     }
     return -w_im_y100(100/(1-x), -x);
@@ -1840,8 +1900,8 @@
 static double relerr(double a, double b) {
   if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) {
     if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) ||
-	(isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) ||
-	(isinf(a) && isinf(b) && a*b < 0))
+        (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) ||
+        (isinf(a) && isinf(b) && a*b < 0))
       return Inf; // "infinite" error
     return 0; // matching infinity/nan results counted as zero error
   }
@@ -1857,173 +1917,173 @@
     printf("############# w(z) tests #############\n");
 #define NTST 57 // define instead of const for C compatibility
     cmplx z[NTST] = {
-      CMPLX(624.2,-0.26123),
-      CMPLX(-0.4,3.),
-      CMPLX(0.6,2.),
-      CMPLX(-1.,1.),
-      CMPLX(-1.,-9.),
-      CMPLX(-1.,9.),
-      CMPLX(-0.0000000234545,1.1234),
-      CMPLX(-3.,5.1),
-      CMPLX(-53,30.1),
-      CMPLX(0.0,0.12345),
-      CMPLX(11,1),
-      CMPLX(-22,-2),
-      CMPLX(9,-28),
-      CMPLX(21,-33),
-      CMPLX(1e5,1e5),
-      CMPLX(1e14,1e14),
-      CMPLX(-3001,-1000),
-      CMPLX(1e160,-1e159),
-      CMPLX(-6.01,0.01),
-      CMPLX(-0.7,-0.7),
-      CMPLX(2.611780000000000e+01, 4.540909610972489e+03),
-      CMPLX(0.8e7,0.3e7),
-      CMPLX(-20,-19.8081),
-      CMPLX(1e-16,-1.1e-16),
-      CMPLX(2.3e-8,1.3e-8),
-      CMPLX(6.3,-1e-13),
-      CMPLX(6.3,1e-20),
-      CMPLX(1e-20,6.3),
-      CMPLX(1e-20,16.3),
-      CMPLX(9,1e-300),
-      CMPLX(6.01,0.11),
-      CMPLX(8.01,1.01e-10),
-      CMPLX(28.01,1e-300),
-      CMPLX(10.01,1e-200),
-      CMPLX(10.01,-1e-200),
-      CMPLX(10.01,0.99e-10),
-      CMPLX(10.01,-0.99e-10),
-      CMPLX(1e-20,7.01),
-      CMPLX(-1,7.01),
-      CMPLX(5.99,7.01),
-      CMPLX(1,0),
-      CMPLX(55,0),
-      CMPLX(-0.1,0),
-      CMPLX(1e-20,0),
-      CMPLX(0,5e-14),
-      CMPLX(0,51),
-      CMPLX(Inf,0),
-      CMPLX(-Inf,0),
-      CMPLX(0,Inf),
-      CMPLX(0,-Inf),
-      CMPLX(Inf,Inf),
-      CMPLX(Inf,-Inf),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(0,NaN),
-      CMPLX(NaN,Inf),
-      CMPLX(Inf,NaN)
+      C(624.2,-0.26123),
+      C(-0.4,3.),
+      C(0.6,2.),
+      C(-1.,1.),
+      C(-1.,-9.),
+      C(-1.,9.),
+      C(-0.0000000234545,1.1234),
+      C(-3.,5.1),
+      C(-53,30.1),
+      C(0.0,0.12345),
+      C(11,1),
+      C(-22,-2),
+      C(9,-28),
+      C(21,-33),
+      C(1e5,1e5),
+      C(1e14,1e14),
+      C(-3001,-1000),
+      C(1e160,-1e159),
+      C(-6.01,0.01),
+      C(-0.7,-0.7),
+      C(2.611780000000000e+01, 4.540909610972489e+03),
+      C(0.8e7,0.3e7),
+      C(-20,-19.8081),
+      C(1e-16,-1.1e-16),
+      C(2.3e-8,1.3e-8),
+      C(6.3,-1e-13),
+      C(6.3,1e-20),
+      C(1e-20,6.3),
+      C(1e-20,16.3),
+      C(9,1e-300),
+      C(6.01,0.11),
+      C(8.01,1.01e-10),
+      C(28.01,1e-300),
+      C(10.01,1e-200),
+      C(10.01,-1e-200),
+      C(10.01,0.99e-10),
+      C(10.01,-0.99e-10),
+      C(1e-20,7.01),
+      C(-1,7.01),
+      C(5.99,7.01),
+      C(1,0),
+      C(55,0),
+      C(-0.1,0),
+      C(1e-20,0),
+      C(0,5e-14),
+      C(0,51),
+      C(Inf,0),
+      C(-Inf,0),
+      C(0,Inf),
+      C(0,-Inf),
+      C(Inf,Inf),
+      C(Inf,-Inf),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(0,NaN),
+      C(NaN,Inf),
+      C(Inf,NaN)
     };
     cmplx w[NTST] = { /* w(z), computed with WolframAlpha
-				   ... note that WolframAlpha is problematic
-				   some of the above inputs, so I had to
-				   use the continued-fraction expansion
-				   in WolframAlpha in some cases, or switch
-				   to Maple */
-      CMPLX(-3.78270245518980507452677445620103199303131110e-7,
-	    0.000903861276433172057331093754199933411710053155),
-      CMPLX(0.1764906227004816847297495349730234591778719532788,
-	    -0.02146550539468457616788719893991501311573031095617),
-      CMPLX(0.2410250715772692146133539023007113781272362309451,
-	    0.06087579663428089745895459735240964093522265589350),
-      CMPLX(0.30474420525691259245713884106959496013413834051768,
-	    -0.20821893820283162728743734725471561394145872072738),
-      CMPLX(7.317131068972378096865595229600561710140617977e34,
-	    8.321873499714402777186848353320412813066170427e34),
-      CMPLX(0.0615698507236323685519612934241429530190806818395,
-	    -0.00676005783716575013073036218018565206070072304635),
-      CMPLX(0.3960793007699874918961319170187598400134746631,
-	    -5.593152259116644920546186222529802777409274656e-9),
-      CMPLX(0.08217199226739447943295069917990417630675021771804,
-	    -0.04701291087643609891018366143118110965272615832184),
-      CMPLX(0.00457246000350281640952328010227885008541748668738,
-	    -0.00804900791411691821818731763401840373998654987934),
-      CMPLX(0.8746342859608052666092782112565360755791467973338452,
-	    0.),
-      CMPLX(0.00468190164965444174367477874864366058339647648741,
-	    0.0510735563901306197993676329845149741675029197050),
-      CMPLX(-0.0023193175200187620902125853834909543869428763219,
-	    -0.025460054739731556004902057663500272721780776336),
-      CMPLX(9.11463368405637174660562096516414499772662584e304,
-	    3.97101807145263333769664875189354358563218932e305),
-      CMPLX(-4.4927207857715598976165541011143706155432296e281,
-	    -2.8019591213423077494444700357168707775769028e281),
-      CMPLX(2.820947917809305132678577516325951485807107151e-6,
-	    2.820947917668257736791638444590253942253354058e-6),
-      CMPLX(2.82094791773878143474039725787438662716372268e-15,
-	    2.82094791773878143474039725773333923127678361e-15),
-      CMPLX(-0.0000563851289696244350147899376081488003110150498,
-	    -0.000169211755126812174631861529808288295454992688),
-      CMPLX(-5.586035480670854326218608431294778077663867e-162,
-	    5.586035480670854326218608431294778077663867e-161),
-      CMPLX(0.00016318325137140451888255634399123461580248456,
-	    -0.095232456573009287370728788146686162555021209999),
-      CMPLX(0.69504753678406939989115375989939096800793577783885,
-	    -1.8916411171103639136680830887017670616339912024317),
-      CMPLX(0.0001242418269653279656612334210746733213167234822,
-	    7.145975826320186888508563111992099992116786763e-7),
-      CMPLX(2.318587329648353318615800865959225429377529825e-8,
-	    6.182899545728857485721417893323317843200933380e-8),
-      CMPLX(-0.0133426877243506022053521927604277115767311800303,
-	    -0.0148087097143220769493341484176979826888871576145),
-      CMPLX(1.00000000000000012412170838050638522857747934,
-	    1.12837916709551279389615890312156495593616433e-16),
-      CMPLX(0.9999999853310704677583504063775310832036830015,
-	    2.595272024519678881897196435157270184030360773e-8),
-      CMPLX(-1.4731421795638279504242963027196663601154624e-15,
-	    0.090727659684127365236479098488823462473074709),
-      CMPLX(5.79246077884410284575834156425396800754409308e-18,
-	    0.0907276596841273652364790985059772809093822374),
-      CMPLX(0.0884658993528521953466533278764830881245144368,
-	    1.37088352495749125283269718778582613192166760e-22),
-      CMPLX(0.0345480845419190424370085249304184266813447878,
-	    2.11161102895179044968099038990446187626075258e-23),
-      CMPLX(6.63967719958073440070225527042829242391918213e-36,
-	    0.0630820900592582863713653132559743161572639353),
-      CMPLX(0.00179435233208702644891092397579091030658500743634,
-	    0.0951983814805270647939647438459699953990788064762),
-      CMPLX(9.09760377102097999924241322094863528771095448e-13,
-	    0.0709979210725138550986782242355007611074966717),
-      CMPLX(7.2049510279742166460047102593255688682910274423e-304,
-	    0.0201552956479526953866611812593266285000876784321),
-      CMPLX(3.04543604652250734193622967873276113872279682e-44,
-	    0.0566481651760675042930042117726713294607499165),
-      CMPLX(3.04543604652250734193622967873276113872279682e-44,
-	    0.0566481651760675042930042117726713294607499165),
-      CMPLX(0.5659928732065273429286988428080855057102069081e-12,
-	    0.056648165176067504292998527162143030538756683302),
-      CMPLX(-0.56599287320652734292869884280802459698927645e-12,
-	    0.0566481651760675042929985271621430305387566833029),
-      CMPLX(0.0796884251721652215687859778119964009569455462,
-	    1.11474461817561675017794941973556302717225126e-22),
-      CMPLX(0.07817195821247357458545539935996687005781943386550,
-	    -0.01093913670103576690766705513142246633056714279654),
-      CMPLX(0.04670032980990449912809326141164730850466208439937,
-	    0.03944038961933534137558064191650437353429669886545),
-      CMPLX(0.36787944117144232159552377016146086744581113103176,
-	    0.60715770584139372911503823580074492116122092866515),
-      CMPLX(0,
-	    0.010259688805536830986089913987516716056946786526145),
-      CMPLX(0.99004983374916805357390597718003655777207908125383,
-	    -0.11208866436449538036721343053869621153527769495574),
-      CMPLX(0.99999999999999999999999999999999999999990000,
-	    1.12837916709551257389615890312154517168802603e-20),
-      CMPLX(0.999999999999943581041645226871305192054749891144158,
-	    0),
-      CMPLX(0.0110604154853277201542582159216317923453996211744250,
-	    0),
-      CMPLX(0,0),
-      CMPLX(0,0),
-      CMPLX(0,0),
-      CMPLX(Inf,0),
-      CMPLX(0,0),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN)
+                                   ... note that WolframAlpha is problematic
+                                   some of the above inputs, so I had to
+                                   use the continued-fraction expansion
+                                   in WolframAlpha in some cases, or switch
+                                   to Maple */
+      C(-3.78270245518980507452677445620103199303131110e-7,
+        0.000903861276433172057331093754199933411710053155),
+      C(0.1764906227004816847297495349730234591778719532788,
+        -0.02146550539468457616788719893991501311573031095617),
+      C(0.2410250715772692146133539023007113781272362309451,
+        0.06087579663428089745895459735240964093522265589350),
+      C(0.30474420525691259245713884106959496013413834051768,
+        -0.20821893820283162728743734725471561394145872072738),
+      C(7.317131068972378096865595229600561710140617977e34,
+        8.321873499714402777186848353320412813066170427e34),
+      C(0.0615698507236323685519612934241429530190806818395,
+        -0.00676005783716575013073036218018565206070072304635),
+      C(0.3960793007699874918961319170187598400134746631,
+        -5.593152259116644920546186222529802777409274656e-9),
+      C(0.08217199226739447943295069917990417630675021771804,
+        -0.04701291087643609891018366143118110965272615832184),
+      C(0.00457246000350281640952328010227885008541748668738,
+        -0.00804900791411691821818731763401840373998654987934),
+      C(0.8746342859608052666092782112565360755791467973338452,
+        0.),
+      C(0.00468190164965444174367477874864366058339647648741,
+        0.0510735563901306197993676329845149741675029197050),
+      C(-0.0023193175200187620902125853834909543869428763219,
+        -0.025460054739731556004902057663500272721780776336),
+      C(9.11463368405637174660562096516414499772662584e304,
+        3.97101807145263333769664875189354358563218932e305),
+      C(-4.4927207857715598976165541011143706155432296e281,
+        -2.8019591213423077494444700357168707775769028e281),
+      C(2.820947917809305132678577516325951485807107151e-6,
+        2.820947917668257736791638444590253942253354058e-6),
+      C(2.82094791773878143474039725787438662716372268e-15,
+        2.82094791773878143474039725773333923127678361e-15),
+      C(-0.0000563851289696244350147899376081488003110150498,
+        -0.000169211755126812174631861529808288295454992688),
+      C(-5.586035480670854326218608431294778077663867e-162,
+        5.586035480670854326218608431294778077663867e-161),
+      C(0.00016318325137140451888255634399123461580248456,
+        -0.095232456573009287370728788146686162555021209999),
+      C(0.69504753678406939989115375989939096800793577783885,
+        -1.8916411171103639136680830887017670616339912024317),
+      C(0.0001242418269653279656612334210746733213167234822,
+        7.145975826320186888508563111992099992116786763e-7),
+      C(2.318587329648353318615800865959225429377529825e-8,
+        6.182899545728857485721417893323317843200933380e-8),
+      C(-0.0133426877243506022053521927604277115767311800303,
+        -0.0148087097143220769493341484176979826888871576145),
+      C(1.00000000000000012412170838050638522857747934,
+        1.12837916709551279389615890312156495593616433e-16),
+      C(0.9999999853310704677583504063775310832036830015,
+        2.595272024519678881897196435157270184030360773e-8),
+      C(-1.4731421795638279504242963027196663601154624e-15,
+        0.090727659684127365236479098488823462473074709),
+      C(5.79246077884410284575834156425396800754409308e-18,
+        0.0907276596841273652364790985059772809093822374),
+      C(0.0884658993528521953466533278764830881245144368,
+        1.37088352495749125283269718778582613192166760e-22),
+      C(0.0345480845419190424370085249304184266813447878,
+        2.11161102895179044968099038990446187626075258e-23),
+      C(6.63967719958073440070225527042829242391918213e-36,
+        0.0630820900592582863713653132559743161572639353),
+      C(0.00179435233208702644891092397579091030658500743634,
+        0.0951983814805270647939647438459699953990788064762),
+      C(9.09760377102097999924241322094863528771095448e-13,
+        0.0709979210725138550986782242355007611074966717),
+      C(7.2049510279742166460047102593255688682910274423e-304,
+        0.0201552956479526953866611812593266285000876784321),
+      C(3.04543604652250734193622967873276113872279682e-44,
+        0.0566481651760675042930042117726713294607499165),
+      C(3.04543604652250734193622967873276113872279682e-44,
+        0.0566481651760675042930042117726713294607499165),
+      C(0.5659928732065273429286988428080855057102069081e-12,
+        0.056648165176067504292998527162143030538756683302),
+      C(-0.56599287320652734292869884280802459698927645e-12,
+        0.0566481651760675042929985271621430305387566833029),
+      C(0.0796884251721652215687859778119964009569455462,
+        1.11474461817561675017794941973556302717225126e-22),
+      C(0.07817195821247357458545539935996687005781943386550,
+        -0.01093913670103576690766705513142246633056714279654),
+      C(0.04670032980990449912809326141164730850466208439937,
+        0.03944038961933534137558064191650437353429669886545),
+      C(0.36787944117144232159552377016146086744581113103176,
+        0.60715770584139372911503823580074492116122092866515),
+      C(0,
+        0.010259688805536830986089913987516716056946786526145),
+      C(0.99004983374916805357390597718003655777207908125383,
+        -0.11208866436449538036721343053869621153527769495574),
+      C(0.99999999999999999999999999999999999999990000,
+        1.12837916709551257389615890312154517168802603e-20),
+      C(0.999999999999943581041645226871305192054749891144158,
+        0),
+      C(0.0110604154853277201542582159216317923453996211744250,
+        0),
+      C(0,0),
+      C(0,0),
+      C(0,0),
+      C(Inf,0),
+      C(0,0),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(NaN,NaN),
+      C(NaN,NaN)
     };
     double errmax = 0;
     for (int i = 0; i < NTST; ++i) {
@@ -2031,8 +2091,8 @@
       double re_err = relerr(creal(w[i]), creal(fw));
       double im_err = relerr(cimag(w[i]), cimag(fw));
       printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n",
-	     creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]),
-	     re_err, im_err);
+             creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]),
+             re_err, im_err);
       if (re_err > errmax) errmax = re_err;
       if (im_err > errmax) errmax = im_err;
     }
@@ -2045,140 +2105,163 @@
   }
   {
 #undef NTST
-#define NTST 33 // define instead of const for C compatibility
+#define NTST 41 // define instead of const for C compatibility
     cmplx z[NTST] = {
-      CMPLX(1,2),
-      CMPLX(-1,2),
-      CMPLX(1,-2),
-      CMPLX(-1,-2),
-      CMPLX(9,-28),
-      CMPLX(21,-33),
-      CMPLX(1e3,1e3),
-      CMPLX(-3001,-1000),
-      CMPLX(1e160,-1e159),
-      CMPLX(5.1e-3, 1e-8),
-      CMPLX(-4.9e-3, 4.95e-3),
-      CMPLX(4.9e-3, 0.5),
-      CMPLX(4.9e-4, -0.5e1),
-      CMPLX(-4.9e-5, -0.5e2),
-      CMPLX(5.1e-3, 0.5),
-      CMPLX(5.1e-4, -0.5e1),
-      CMPLX(-5.1e-5, -0.5e2),
-      CMPLX(1e-6,2e-6),
-      CMPLX(0,2e-6),
-      CMPLX(0,2),
-      CMPLX(0,20),
-      CMPLX(0,200),
-      CMPLX(Inf,0),
-      CMPLX(-Inf,0),
-      CMPLX(0,Inf),
-      CMPLX(0,-Inf),
-      CMPLX(Inf,Inf),
-      CMPLX(Inf,-Inf),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(0,NaN),
-      CMPLX(NaN,Inf),
-      CMPLX(Inf,NaN)
+      C(1,2),
+      C(-1,2),
+      C(1,-2),
+      C(-1,-2),
+      C(9,-28),
+      C(21,-33),
+      C(1e3,1e3),
+      C(-3001,-1000),
+      C(1e160,-1e159),
+      C(5.1e-3, 1e-8),
+      C(-4.9e-3, 4.95e-3),
+      C(4.9e-3, 0.5),
+      C(4.9e-4, -0.5e1),
+      C(-4.9e-5, -0.5e2),
+      C(5.1e-3, 0.5),
+      C(5.1e-4, -0.5e1),
+      C(-5.1e-5, -0.5e2),
+      C(1e-6,2e-6),
+      C(0,2e-6),
+      C(0,2),
+      C(0,20),
+      C(0,200),
+      C(Inf,0),
+      C(-Inf,0),
+      C(0,Inf),
+      C(0,-Inf),
+      C(Inf,Inf),
+      C(Inf,-Inf),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(0,NaN),
+      C(NaN,Inf),
+      C(Inf,NaN),
+      C(1e-3,NaN),
+      C(7e-2,7e-2),
+      C(7e-2,-7e-4),
+      C(-9e-2,7e-4),
+      C(-9e-2,9e-2),
+      C(-7e-4,9e-2),
+      C(7e-2,0.9e-2),
+      C(7e-2,1.1e-2)
     };
     cmplx w[NTST] = { // erf(z[i]), evaluated with Maple
-      CMPLX(-0.5366435657785650339917955593141927494421,
-	    -5.049143703447034669543036958614140565553),
-      CMPLX(0.5366435657785650339917955593141927494421,
-	    -5.049143703447034669543036958614140565553),
-      CMPLX(-0.5366435657785650339917955593141927494421,
-	    5.049143703447034669543036958614140565553),
-      CMPLX(0.5366435657785650339917955593141927494421,
-	    5.049143703447034669543036958614140565553),
-      CMPLX(0.3359473673830576996788000505817956637777e304,
-	    -0.1999896139679880888755589794455069208455e304),
-      CMPLX(0.3584459971462946066523939204836760283645e278,
-	    0.3818954885257184373734213077678011282505e280),
-      CMPLX(0.9996020422657148639102150147542224526887,
-	    0.00002801044116908227889681753993542916894856),
-      CMPLX(-1, 0),
-      CMPLX(1, 0),
-      CMPLX(0.005754683859034800134412990541076554934877,
-	    0.1128349818335058741511924929801267822634e-7),
-      CMPLX(-0.005529149142341821193633460286828381876955,
-	    0.005585388387864706679609092447916333443570),
-      CMPLX(0.007099365669981359632319829148438283865814,
-	    0.6149347012854211635026981277569074001219),
-      CMPLX(0.3981176338702323417718189922039863062440e8,
-	    -0.8298176341665249121085423917575122140650e10),
-      CMPLX(-Inf,
-	    -Inf),
-      CMPLX(0.007389128308257135427153919483147229573895,
-	    0.6149332524601658796226417164791221815139),
-      CMPLX(0.4143671923267934479245651547534414976991e8,
-	    -0.8298168216818314211557046346850921446950e10),
-      CMPLX(-Inf,
-	    -Inf),
-      CMPLX(0.1128379167099649964175513742247082845155e-5,
-	    0.2256758334191777400570377193451519478895e-5),
-      CMPLX(0,
-	    0.2256758334194034158904576117253481476197e-5),
-      CMPLX(0,
-	    18.56480241457555259870429191324101719886),
-      CMPLX(0,
-	    0.1474797539628786202447733153131835124599e173),
-      CMPLX(0,
-	    Inf),
-      CMPLX(1,0),
-      CMPLX(-1,0),
-      CMPLX(0,Inf),
-      CMPLX(0,-Inf),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(0,NaN),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,NaN)
+      C(-0.5366435657785650339917955593141927494421,
+        -5.049143703447034669543036958614140565553),
+      C(0.5366435657785650339917955593141927494421,
+        -5.049143703447034669543036958614140565553),
+      C(-0.5366435657785650339917955593141927494421,
+        5.049143703447034669543036958614140565553),
+      C(0.5366435657785650339917955593141927494421,
+        5.049143703447034669543036958614140565553),
+      C(0.3359473673830576996788000505817956637777e304,
+        -0.1999896139679880888755589794455069208455e304),
+      C(0.3584459971462946066523939204836760283645e278,
+        0.3818954885257184373734213077678011282505e280),
+      C(0.9996020422657148639102150147542224526887,
+        0.00002801044116908227889681753993542916894856),
+      C(-1, 0),
+      C(1, 0),
+      C(0.005754683859034800134412990541076554934877,
+        0.1128349818335058741511924929801267822634e-7),
+      C(-0.005529149142341821193633460286828381876955,
+        0.005585388387864706679609092447916333443570),
+      C(0.007099365669981359632319829148438283865814,
+        0.6149347012854211635026981277569074001219),
+      C(0.3981176338702323417718189922039863062440e8,
+        -0.8298176341665249121085423917575122140650e10),
+      C(-Inf,
+        -Inf),
+      C(0.007389128308257135427153919483147229573895,
+        0.6149332524601658796226417164791221815139),
+      C(0.4143671923267934479245651547534414976991e8,
+        -0.8298168216818314211557046346850921446950e10),
+      C(-Inf,
+        -Inf),
+      C(0.1128379167099649964175513742247082845155e-5,
+        0.2256758334191777400570377193451519478895e-5),
+      C(0,
+        0.2256758334194034158904576117253481476197e-5),
+      C(0,
+        18.56480241457555259870429191324101719886),
+      C(0,
+        0.1474797539628786202447733153131835124599e173),
+      C(0,
+        Inf),
+      C(1,0),
+      C(-1,0),
+      C(0,Inf),
+      C(0,-Inf),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(0,NaN),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(NaN,NaN),
+      C(0.07924380404615782687930591956705225541145,
+        0.07872776218046681145537914954027729115247),
+      C(0.07885775828512276968931773651224684454495,
+        -0.0007860046704118224342390725280161272277506),
+      C(-0.1012806432747198859687963080684978759881,
+        0.0007834934747022035607566216654982820299469),
+      C(-0.1020998418798097910247132140051062512527,
+        0.1010030778892310851309082083238896270340),
+      C(-0.0007962891763147907785684591823889484764272,
+        0.1018289385936278171741809237435404896152),
+      C(0.07886408666470478681566329888615410479530,
+        0.01010604288780868961492224347707949372245),
+      C(0.07886723099940260286824654364807981336591,
+        0.01235199327873258197931147306290916629654)
     };
-#define TST(f,isc)							\
-    printf("############# " #f "(z) tests #############\n");		\
-    double errmax = 0;							\
-    for (int i = 0; i < NTST; ++i) {					\
-      cmplx fw = FADDEEVA(f)(z[i],0.);			\
-      double re_err = relerr(creal(w[i]), creal(fw));			\
-      double im_err = relerr(cimag(w[i]), cimag(fw));			\
+#define TST(f,isc)                                                      \
+    printf("############# " #f "(z) tests #############\n");            \
+    double errmax = 0;                                                  \
+    for (int i = 0; i < NTST; ++i) {                                    \
+      cmplx fw = FADDEEVA(f)(z[i],0.);                  \
+      double re_err = relerr(creal(w[i]), creal(fw));                   \
+      double im_err = relerr(cimag(w[i]), cimag(fw));                   \
       printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \
-	     creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \
-	     re_err, im_err);						\
-      if (re_err > errmax) errmax = re_err;				\
-      if (im_err > errmax) errmax = im_err;				\
-    }									\
-    if (errmax > 1e-13) {						\
-      printf("FAILURE -- relative error %g too large!\n", errmax);	\
-      return 1;								\
-    }									\
-    printf("Checking " #f "(x) special case...\n");			\
-    for (int i = 0; i < 10000; ++i) {					\
-      double x = pow(10., -300. + i * 600. / (10000 - 1));		\
-      double re_err = relerr(FADDEEVA_RE(f)(x),			        \
-			     creal(FADDEEVA(f)(CMPLX(x,x*isc),0.)));	\
-      if (re_err > errmax) errmax = re_err;				\
-      re_err = relerr(FADDEEVA_RE(f)(-x),				\
-		      creal(FADDEEVA(f)(CMPLX(-x,x*isc),0.)));		\
-      if (re_err > errmax) errmax = re_err;				\
-    }									\
-    {									\
-      double re_err = relerr(FADDEEVA_RE(f)(Inf),			\
-			     creal(FADDEEVA(f)(CMPLX(Inf,0.),0.)));	\
-      if (re_err > errmax) errmax = re_err;				\
-      re_err = relerr(FADDEEVA_RE(f)(-Inf),				\
-		      creal(FADDEEVA(f)(CMPLX(-Inf,0.),0.)));		\
-      if (re_err > errmax) errmax = re_err;				\
-      re_err = relerr(FADDEEVA_RE(f)(NaN),				\
-		      creal(FADDEEVA(f)(CMPLX(NaN,0.),0.)));		\
-      if (re_err > errmax) errmax = re_err;				\
-    }									\
-    if (errmax > 1e-13) {						\
-      printf("FAILURE -- relative error %g too large!\n", errmax);	\
-      return 1;								\
-    }									\
-    printf("SUCCESS (max relative error = %g)\n", errmax);		\
+             creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \
+             re_err, im_err);                                           \
+      if (re_err > errmax) errmax = re_err;                             \
+      if (im_err > errmax) errmax = im_err;                             \
+    }                                                                   \
+    if (errmax > 1e-13) {                                               \
+      printf("FAILURE -- relative error %g too large!\n", errmax);      \
+      return 1;                                                         \
+    }                                                                   \
+    printf("Checking " #f "(x) special case...\n");                     \
+    for (int i = 0; i < 10000; ++i) {                                   \
+      double x = pow(10., -300. + i * 600. / (10000 - 1));              \
+      double re_err = relerr(FADDEEVA_RE(f)(x),                         \
+                             creal(FADDEEVA(f)(C(x,x*isc),0.)));        \
+      if (re_err > errmax) errmax = re_err;                             \
+      re_err = relerr(FADDEEVA_RE(f)(-x),                               \
+                      creal(FADDEEVA(f)(C(-x,x*isc),0.)));              \
+      if (re_err > errmax) errmax = re_err;                             \
+    }                                                                   \
+    {                                                                   \
+      double re_err = relerr(FADDEEVA_RE(f)(Inf),                       \
+                             creal(FADDEEVA(f)(C(Inf,0.),0.))); \
+      if (re_err > errmax) errmax = re_err;                             \
+      re_err = relerr(FADDEEVA_RE(f)(-Inf),                             \
+                      creal(FADDEEVA(f)(C(-Inf,0.),0.)));               \
+      if (re_err > errmax) errmax = re_err;                             \
+      re_err = relerr(FADDEEVA_RE(f)(NaN),                              \
+                      creal(FADDEEVA(f)(C(NaN,0.),0.)));                \
+      if (re_err > errmax) errmax = re_err;                             \
+    }                                                                   \
+    if (errmax > 1e-13) {                                               \
+      printf("FAILURE -- relative error %g too large!\n", errmax);      \
+      return 1;                                                         \
+    }                                                                   \
+    printf("SUCCESS (max relative error = %g)\n", errmax);              \
     if (errmax > errmax_all) errmax_all = errmax
 
     TST(erf, 1e-20);
@@ -2188,10 +2271,10 @@
     // be sufficient to make sure I didn't screw up the signs or something
 #undef NTST
 #define NTST 1 // define instead of const for C compatibility
-    cmplx z[NTST] = { CMPLX(1.234,0.5678) };
+    cmplx z[NTST] = { C(1.234,0.5678) };
     cmplx w[NTST] = { // erfi(z[i]), computed with Maple
-      CMPLX(1.081032284405373149432716643834106923212,
-		      1.926775520840916645838949402886591180834)
+      C(1.081032284405373149432716643834106923212,
+        1.926775520840916645838949402886591180834)
     };
     TST(erfi, 0);
   }
@@ -2200,10 +2283,10 @@
     // be sufficient to make sure I didn't screw up the signs or something
 #undef NTST
 #define NTST 1 // define instead of const for C compatibility
-    cmplx z[NTST] = { CMPLX(1.234,0.5678) };
+    cmplx z[NTST] = { C(1.234,0.5678) };
     cmplx w[NTST] = { // erfcx(z[i]), computed with Maple
-      CMPLX(0.3382187479799972294747793561190487832579,
-		      -0.1116077470811648467464927471872945833154)
+      C(0.3382187479799972294747793561190487832579,
+        -0.1116077470811648467464927471872945833154)
     };
     TST(erfcx, 0);
   }
@@ -2211,214 +2294,217 @@
 #undef NTST
 #define NTST 30 // define instead of const for C compatibility
     cmplx z[NTST] = {
-      CMPLX(1,2),
-      CMPLX(-1,2),
-      CMPLX(1,-2),
-      CMPLX(-1,-2),
-      CMPLX(9,-28),
-      CMPLX(21,-33),
-      CMPLX(1e3,1e3),
-      CMPLX(-3001,-1000),
-      CMPLX(1e160,-1e159),
-      CMPLX(5.1e-3, 1e-8),
-      CMPLX(0,2e-6),
-      CMPLX(0,2),
-      CMPLX(0,20),
-      CMPLX(0,200),
-      CMPLX(2e-6,0),
-      CMPLX(2,0),
-      CMPLX(20,0),
-      CMPLX(200,0),
-      CMPLX(Inf,0),
-      CMPLX(-Inf,0),
-      CMPLX(0,Inf),
-      CMPLX(0,-Inf),
-      CMPLX(Inf,Inf),
-      CMPLX(Inf,-Inf),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(0,NaN),
-      CMPLX(NaN,Inf),
-      CMPLX(Inf,NaN),
-      CMPLX(88,0)
+      C(1,2),
+      C(-1,2),
+      C(1,-2),
+      C(-1,-2),
+      C(9,-28),
+      C(21,-33),
+      C(1e3,1e3),
+      C(-3001,-1000),
+      C(1e160,-1e159),
+      C(5.1e-3, 1e-8),
+      C(0,2e-6),
+      C(0,2),
+      C(0,20),
+      C(0,200),
+      C(2e-6,0),
+      C(2,0),
+      C(20,0),
+      C(200,0),
+      C(Inf,0),
+      C(-Inf,0),
+      C(0,Inf),
+      C(0,-Inf),
+      C(Inf,Inf),
+      C(Inf,-Inf),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(0,NaN),
+      C(NaN,Inf),
+      C(Inf,NaN),
+      C(88,0)
     };
     cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple
-      CMPLX(1.536643565778565033991795559314192749442,
-	    5.049143703447034669543036958614140565553),
-      CMPLX(0.4633564342214349660082044406858072505579,
-	    5.049143703447034669543036958614140565553),
-      CMPLX(1.536643565778565033991795559314192749442,
-	    -5.049143703447034669543036958614140565553),
-      CMPLX(0.4633564342214349660082044406858072505579,
-	    -5.049143703447034669543036958614140565553),
-      CMPLX(-0.3359473673830576996788000505817956637777e304,
-	    0.1999896139679880888755589794455069208455e304),
-      CMPLX(-0.3584459971462946066523939204836760283645e278,
-	    -0.3818954885257184373734213077678011282505e280),
-      CMPLX(0.0003979577342851360897849852457775473112748,
-	    -0.00002801044116908227889681753993542916894856),
-      CMPLX(2, 0),
-      CMPLX(0, 0),
-      CMPLX(0.9942453161409651998655870094589234450651,
-	    -0.1128349818335058741511924929801267822634e-7),
-      CMPLX(1,
-	    -0.2256758334194034158904576117253481476197e-5),
-      CMPLX(1,
-	    -18.56480241457555259870429191324101719886),
-      CMPLX(1,
-	    -0.1474797539628786202447733153131835124599e173),
-      CMPLX(1, -Inf),
-      CMPLX(0.9999977432416658119838633199332831406314,
-	    0),
-      CMPLX(0.004677734981047265837930743632747071389108,
-	    0),
-      CMPLX(0.5395865611607900928934999167905345604088e-175,
-	    0),
-      CMPLX(0, 0),
-      CMPLX(0, 0),
-      CMPLX(2, 0),
-      CMPLX(1, -Inf),
-      CMPLX(1, Inf),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, 0),
-      CMPLX(1, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(0,0)
+      C(1.536643565778565033991795559314192749442,
+        5.049143703447034669543036958614140565553),
+      C(0.4633564342214349660082044406858072505579,
+        5.049143703447034669543036958614140565553),
+      C(1.536643565778565033991795559314192749442,
+        -5.049143703447034669543036958614140565553),
+      C(0.4633564342214349660082044406858072505579,
+        -5.049143703447034669543036958614140565553),
+      C(-0.3359473673830576996788000505817956637777e304,
+        0.1999896139679880888755589794455069208455e304),
+      C(-0.3584459971462946066523939204836760283645e278,
+        -0.3818954885257184373734213077678011282505e280),
+      C(0.0003979577342851360897849852457775473112748,
+        -0.00002801044116908227889681753993542916894856),
+      C(2, 0),
+      C(0, 0),
+      C(0.9942453161409651998655870094589234450651,
+        -0.1128349818335058741511924929801267822634e-7),
+      C(1,
+        -0.2256758334194034158904576117253481476197e-5),
+      C(1,
+        -18.56480241457555259870429191324101719886),
+      C(1,
+        -0.1474797539628786202447733153131835124599e173),
+      C(1, -Inf),
+      C(0.9999977432416658119838633199332831406314,
+        0),
+      C(0.004677734981047265837930743632747071389108,
+        0),
+      C(0.5395865611607900928934999167905345604088e-175,
+        0),
+      C(0, 0),
+      C(0, 0),
+      C(2, 0),
+      C(1, -Inf),
+      C(1, Inf),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(NaN, 0),
+      C(1, NaN),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(0,0)
     };
     TST(erfc, 1e-20);
   }
   {
 #undef NTST
-#define NTST 47 // define instead of const for C compatibility
+#define NTST 48 // define instead of const for C compatibility
     cmplx z[NTST] = {
-      CMPLX(2,1),
-      CMPLX(-2,1),
-      CMPLX(2,-1),
-      CMPLX(-2,-1),
-      CMPLX(-28,9),
-      CMPLX(33,-21),
-      CMPLX(1e3,1e3),
-      CMPLX(-1000,-3001),
-      CMPLX(1e-8, 5.1e-3),
-      CMPLX(4.95e-3, -4.9e-3),
-      CMPLX(0.5, 4.9e-3),
-      CMPLX(-0.5e1, 4.9e-4),
-      CMPLX(-0.5e2, -4.9e-5),
-      CMPLX(0.5e3, 4.9e-6),
-      CMPLX(0.5, 5.1e-3),
-      CMPLX(-0.5e1, 5.1e-4),
-      CMPLX(-0.5e2, -5.1e-5),
-      CMPLX(1e-6,2e-6),
-      CMPLX(2e-6,0),
-      CMPLX(2,0),
-      CMPLX(20,0),
-      CMPLX(200,0),
-      CMPLX(0,4.9e-3),
-      CMPLX(0,-5.1e-3),
-      CMPLX(0,2e-6),
-      CMPLX(0,-2),
-      CMPLX(0,20),
-      CMPLX(0,-200),
-      CMPLX(Inf,0),
-      CMPLX(-Inf,0),
-      CMPLX(0,Inf),
-      CMPLX(0,-Inf),
-      CMPLX(Inf,Inf),
-      CMPLX(Inf,-Inf),
-      CMPLX(NaN,NaN),
-      CMPLX(NaN,0),
-      CMPLX(0,NaN),
-      CMPLX(NaN,Inf),
-      CMPLX(Inf,NaN),
-      CMPLX(39, 6.4e-5),
-      CMPLX(41, 6.09e-5),
-      CMPLX(4.9e7, 5e-11),
-      CMPLX(5.1e7, 4.8e-11),
-      CMPLX(1e9, 2.4e-12),
-      CMPLX(1e11, 2.4e-14),
-      CMPLX(1e13, 2.4e-16),
-      CMPLX(1e300, 2.4e-303)
+      C(2,1),
+      C(-2,1),
+      C(2,-1),
+      C(-2,-1),
+      C(-28,9),
+      C(33,-21),
+      C(1e3,1e3),
+      C(-1000,-3001),
+      C(1e-8, 5.1e-3),
+      C(4.95e-3, -4.9e-3),
+      C(5.1e-3, 5.1e-3),
+      C(0.5, 4.9e-3),
+      C(-0.5e1, 4.9e-4),
+      C(-0.5e2, -4.9e-5),
+      C(0.5e3, 4.9e-6),
+      C(0.5, 5.1e-3),
+      C(-0.5e1, 5.1e-4),
+      C(-0.5e2, -5.1e-5),
+      C(1e-6,2e-6),
+      C(2e-6,0),
+      C(2,0),
+      C(20,0),
+      C(200,0),
+      C(0,4.9e-3),
+      C(0,-5.1e-3),
+      C(0,2e-6),
+      C(0,-2),
+      C(0,20),
+      C(0,-200),
+      C(Inf,0),
+      C(-Inf,0),
+      C(0,Inf),
+      C(0,-Inf),
+      C(Inf,Inf),
+      C(Inf,-Inf),
+      C(NaN,NaN),
+      C(NaN,0),
+      C(0,NaN),
+      C(NaN,Inf),
+      C(Inf,NaN),
+      C(39, 6.4e-5),
+      C(41, 6.09e-5),
+      C(4.9e7, 5e-11),
+      C(5.1e7, 4.8e-11),
+      C(1e9, 2.4e-12),
+      C(1e11, 2.4e-14),
+      C(1e13, 2.4e-16),
+      C(1e300, 2.4e-303)
     };
     cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple
-      CMPLX(0.1635394094345355614904345232875688576839,
-	    -0.1531245755371229803585918112683241066853),
-      CMPLX(-0.1635394094345355614904345232875688576839,
-	    -0.1531245755371229803585918112683241066853),
-      CMPLX(0.1635394094345355614904345232875688576839,
-	    0.1531245755371229803585918112683241066853),
-      CMPLX(-0.1635394094345355614904345232875688576839,
-	    0.1531245755371229803585918112683241066853),
-      CMPLX(-0.01619082256681596362895875232699626384420,
-	    -0.005210224203359059109181555401330902819419),
-      CMPLX(0.01078377080978103125464543240346760257008,
-	    0.006866888783433775382193630944275682670599),
-      CMPLX(-0.5808616819196736225612296471081337245459,
-	    0.6688593905505562263387760667171706325749),
-      CMPLX(Inf,
-	    -Inf),
-      CMPLX(0.1000052020902036118082966385855563526705e-7,
-	    0.005100088434920073153418834680320146441685),
-      CMPLX(0.004950156837581592745389973960217444687524,
-	    -0.004899838305155226382584756154100963570500),
-      CMPLX(0.4244534840871830045021143490355372016428,
-	    0.002820278933186814021399602648373095266538),
-      CMPLX(-0.1021340733271046543881236523269967674156,
-	    -0.00001045696456072005761498961861088944159916),
-      CMPLX(-0.01000200120119206748855061636187197886859,
-	    0.9805885888237419500266621041508714123763e-8),
-      CMPLX(0.001000002000012000023960527532953151819595,
-	    -0.9800058800588007290937355024646722133204e-11),
-      CMPLX(0.4244549085628511778373438768121222815752,
-	    0.002935393851311701428647152230552122898291),
-      CMPLX(-0.1021340732357117208743299813648493928105,
-	    -0.00001088377943049851799938998805451564893540),
-      CMPLX(-0.01000200120119126652710792390331206563616,
-	    0.1020612612857282306892368985525393707486e-7),
-      CMPLX(0.1000000000007333333333344266666666664457e-5,
-	    0.2000000000001333333333323199999999978819e-5),
-      CMPLX(0.1999999999994666666666675199999999990248e-5,
-	    0),
-      CMPLX(0.3013403889237919660346644392864226952119,
-	    0),
-      CMPLX(0.02503136792640367194699495234782353186858,
-	    0),
-      CMPLX(0.002500031251171948248596912483183760683918,
-	    0),
-      CMPLX(0,0.004900078433419939164774792850907128053308),
-      CMPLX(0,-0.005100088434920074173454208832365950009419),
-      CMPLX(0,0.2000000000005333333333341866666666676419e-5),
-      CMPLX(0,-48.16001211429122974789822893525016528191),
-      CMPLX(0,0.4627407029504443513654142715903005954668e174),
-      CMPLX(0,-Inf),
-      CMPLX(0,0),
-      CMPLX(-0,0),
-      CMPLX(0, Inf),
-      CMPLX(0, -Inf),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, 0),
-      CMPLX(0, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(NaN, NaN),
-      CMPLX(0.01282473148489433743567240624939698290584,
-	    -0.2105957276516618621447832572909153498104e-7),
-      CMPLX(0.01219875253423634378984109995893708152885,
-	    -0.1813040560401824664088425926165834355953e-7),
-      CMPLX(0.1020408163265306334945473399689037886997e-7,
-	    -0.1041232819658476285651490827866174985330e-25),
-      CMPLX(0.9803921568627452865036825956835185367356e-8,
-	    -0.9227220299884665067601095648451913375754e-26),
-      CMPLX(0.5000000000000000002500000000000000003750e-9,
-	    -0.1200000000000000001800000188712838420241e-29),
-      CMPLX(5.00000000000000000000025000000000000000000003e-12,
-	    -1.20000000000000000000018000000000000000000004e-36),
-      CMPLX(5.00000000000000000000000002500000000000000000e-14,
-	    -1.20000000000000000000000001800000000000000000e-42),
-      CMPLX(5e-301, 0)
+      C(0.1635394094345355614904345232875688576839,
+        -0.1531245755371229803585918112683241066853),
+      C(-0.1635394094345355614904345232875688576839,
+        -0.1531245755371229803585918112683241066853),
+      C(0.1635394094345355614904345232875688576839,
+        0.1531245755371229803585918112683241066853),
+      C(-0.1635394094345355614904345232875688576839,
+        0.1531245755371229803585918112683241066853),
+      C(-0.01619082256681596362895875232699626384420,
+        -0.005210224203359059109181555401330902819419),
+      C(0.01078377080978103125464543240346760257008,
+        0.006866888783433775382193630944275682670599),
+      C(-0.5808616819196736225612296471081337245459,
+        0.6688593905505562263387760667171706325749),
+      C(Inf,
+        -Inf),
+      C(0.1000052020902036118082966385855563526705e-7,
+        0.005100088434920073153418834680320146441685),
+      C(0.004950156837581592745389973960217444687524,
+        -0.004899838305155226382584756154100963570500),
+      C(0.005100176864319675957314822982399286703798,
+        0.005099823128319785355949825238269336481254),
+      C(0.4244534840871830045021143490355372016428,
+        0.002820278933186814021399602648373095266538),
+      C(-0.1021340733271046543881236523269967674156,
+        -0.00001045696456072005761498961861088944159916),
+      C(-0.01000200120119206748855061636187197886859,
+        0.9805885888237419500266621041508714123763e-8),
+      C(0.001000002000012000023960527532953151819595,
+        -0.9800058800588007290937355024646722133204e-11),
+      C(0.4244549085628511778373438768121222815752,
+        0.002935393851311701428647152230552122898291),
+      C(-0.1021340732357117208743299813648493928105,
+        -0.00001088377943049851799938998805451564893540),
+      C(-0.01000200120119126652710792390331206563616,
+        0.1020612612857282306892368985525393707486e-7),
+      C(0.1000000000007333333333344266666666664457e-5,
+        0.2000000000001333333333323199999999978819e-5),
+      C(0.1999999999994666666666675199999999990248e-5,
+        0),
+      C(0.3013403889237919660346644392864226952119,
+        0),
+      C(0.02503136792640367194699495234782353186858,
+        0),
+      C(0.002500031251171948248596912483183760683918,
+        0),
+      C(0,0.004900078433419939164774792850907128053308),
+      C(0,-0.005100088434920074173454208832365950009419),
+      C(0,0.2000000000005333333333341866666666676419e-5),
+      C(0,-48.16001211429122974789822893525016528191),
+      C(0,0.4627407029504443513654142715903005954668e174),
+      C(0,-Inf),
+      C(0,0),
+      C(-0,0),
+      C(0, Inf),
+      C(0, -Inf),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(NaN, 0),
+      C(0, NaN),
+      C(NaN, NaN),
+      C(NaN, NaN),
+      C(0.01282473148489433743567240624939698290584,
+        -0.2105957276516618621447832572909153498104e-7),
+      C(0.01219875253423634378984109995893708152885,
+        -0.1813040560401824664088425926165834355953e-7),
+      C(0.1020408163265306334945473399689037886997e-7,
+        -0.1041232819658476285651490827866174985330e-25),
+      C(0.9803921568627452865036825956835185367356e-8,
+        -0.9227220299884665067601095648451913375754e-26),
+      C(0.5000000000000000002500000000000000003750e-9,
+        -0.1200000000000000001800000188712838420241e-29),
+      C(5.00000000000000000000025000000000000000000003e-12,
+        -1.20000000000000000000018000000000000000000004e-36),
+      C(5.00000000000000000000000002500000000000000000e-14,
+        -1.20000000000000000000000001800000000000000000e-42),
+      C(5e-301, 0)
     };
     TST(Dawson, 1e-20);
   }