Faddeeva.cpp

//  -*- 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
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 * 
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
 */

/* (Note that this file can be compiled with either C++, in which
    case it uses C++ std::complex<double>, or C, in which case it
    uses C99 double complex.) */

/* Available at: http://ab-initio.mit.edu/Faddeeva

   Computes various error functions (erf, erfc, erfi, erfcx), 
   including the Dawson integral, in the complex plane, based
   on algorithms for the computation of the Faddeeva function 
              w(z) = exp(-z^2) * erfc(-i*z).
   Given w(z), the error functions are mostly straightforward
   to compute, except for certain regions where we have to
   switch to Taylor expansions to avoid cancellation errors
   [e.g. near the origin for erf(z)].

   To compute the Faddeeva function, we use a combination of two
   algorithms:

   For sufficiently large |z|, we use a continued-fraction expansion
   for w(z) similar to those described in:

      Walter Gautschi, "Efficient computation of the complex error
      function," SIAM J. Numer. Anal. 7(1), pp. 187-198 (1970)

      G. P. M. Poppe and C. M. J. Wijers, "More efficient computation
      of the complex error function," ACM Trans. Math. Soft. 16(1),
      pp. 38-46 (1990).

   Unlike those papers, however, we switch to a completely different
   algorithm for smaller |z|:

      Mofreh R. Zaghloul and Ahmed N. Ali, "Algorithm 916: Computing the
      Faddeyeva and Voigt Functions," ACM Trans. Math. Soft. 38(2), 15
      (2011).

   (I initially used this algorithm for all z, but it turned out to be
    significantly slower than the continued-fraction expansion for
    larger |z|.  On the other hand, it is competitive for smaller |z|, 
    and is significantly more accurate than the Poppe & Wijers code
    in some regions, e.g. in the vicinity of z=1+1i.)

   Note that this is an INDEPENDENT RE-IMPLEMENTATION of these algorithms,
   based on the description in the papers ONLY.  In particular, I did
   not refer to the authors' Fortran or Matlab implementations, respectively,
   (which are under restrictive ACM copyright terms and therefore unusable
    in free/open-source software).

   Steven G. Johnson, Massachusetts Institute of Technology
   http://math.mit.edu/~stevenj
   October 2012.

    -- Note that Algorithm 916 assumes that the erfc(x) function, 
       or rather the scaled function erfcx(x) = exp(x*x)*erfc(x),
       is supplied for REAL arguments x.   I originally used an
       erfcx routine derived from DERFC in SLATEC, but I have
       since replaced it with a much faster routine written by
       me which uses a combination of continued-fraction expansions
       and a lookup table of Chebyshev polynomials.  For speed,
       I implemented a similar algorithm for Im[w(x)] of real x,
       since this comes up frequently in the other error functions.

   A small test program is included the end, which checks
   the w(z) etc. results against several known values.  To compile
   the test function, compile with -DTEST_FADDEEVA (that is,
   #define TEST_FADDEEVA).

   If HAVE_CONFIG_H is #defined (e.g. by compiling with -DHAVE_CONFIG_H),
   then we #include "config.h", which is assumed to be a GNU autoconf-style
   header defining HAVE_* macros to indicate the presence of features. In
   particular, if HAVE_ISNAN and HAVE_ISINF are #defined, we use those
   functions in math.h instead of defining our own, and if HAVE_ERF and/or
   HAVE_ERFC are defined we use those functions from <cmath> for erf and
   erfc of real arguments, respectively, instead of defining our own.

   REVISION HISTORY:
       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)
      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.
      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).
      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.
      29 October 2012: Replace SLATEC-derived erfcx routine with
                       completely rewritten code by me, using a very
                       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.
      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.
      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.
     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.
     18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson)
*/

/////////////////////////////////////////////////////////////////////////
/* If this file is compiled as a part of a larger project,
   support using an autoconf-style config.h header file
   (with various "HAVE_*" #defines to indicate features)
   if HAVE_CONFIG_H is #defined (in GNU autotools style). */

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif

/////////////////////////////////////////////////////////////////////////
// macros to allow us to use either C++ or C (with C99 features)

#ifdef __cplusplus

#  include "Faddeeva.hh"

#  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
#  define cexp(z) exp(z)
#  define creal(z) real(z)
#  define cimag(z) imag(z)
#  define cpolar(r,t) polar(r,t)

#  define C(a,b) cmplx(a,b)

#  define FADDEEVA(name) Faddeeva::name
#  define FADDEEVA_RE(name) Faddeeva::name

// isnan/isinf were introduced in C++11
#  if (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF))
static inline bool my_isnan(double x) { return x != x; }
#    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 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 C(r * cos(t), r * sin(t));
}

#endif // !__cplusplus, i.e. pure C (requires C99 features)

/////////////////////////////////////////////////////////////////////////
// Auxiliary routines to compute other special functions based on w(z)

// compute erfcx(z) = exp(z^2) erfz(z)
cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
{
  return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
}

// compute the error function erf(x)
double FADDEEVA_RE(erf)(double x)
{
#if !defined(__cplusplus)
  return erf(x); // C99 supplies erf in math.h
#elif (__cplusplus >= 201103L) || defined(HAVE_ERF)
  return ::erf(x); // C++11 supplies std::erf in cmath
#else
  double mx2 = -x*x;
  if (mx2 < -750) // underflow
    return (x >= 0 ? 1.0 : -1.0);

  if (x >= 0) {
    if (x < 8e-2) goto taylor;
    return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
  }
  else { // x < 0
    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 - x^6/42 + x^8/216 + ...)
 taylor:
  return x * (1.1283791670955125739
              + mx2 * (0.37612638903183752464
                       + mx2 * (0.11283791670955125739
                                + mx2 * (0.026866170645131251760
                                         + mx2 * 0.0052239776254421878422))));
#endif
}

// compute the error function erf(z)
cmplx FADDEEVA(erf)(cmplx z, double relerr)
{
  double x = creal(z), y = cimag(z);

  if (y == 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 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
    return (x >= 0 ? 1.0 : -1.0);

  /* Handle positive and negative x via different formulas,
     using the mirror symmetries of w, to avoid overflow/underflow
     problems from multiplying exponentially large and small quantities. */
  if (x >= 0) {
    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) *
      (C(cos(mIm_z2), sin(mIm_z2))
       * FADDEEVA(w)(C(-y,x), relerr));
  }
  else { // x < 0
    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 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) *
      (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 - z^6/42 + z^8/216 + ...)
 taylor:
  {
    cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
    return z * (1.1283791670955125739
                + 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 + ... 
              - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
     where:
        erf(iy) = exp(y^2) * Im[w(y)]
  */
 taylor_erfi:
  {
    double x2 = x*x, y2 = y*y;
    double expy2 = exp(y2);
    return C
      (expy2 * x * (1.1283791670955125739
                    - 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))));
  }
}

// erfi(z) = -i erf(iz)
cmplx FADDEEVA(erfi)(cmplx z, double relerr)
{
  cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr);
  return C(cimag(e), -creal(e));
}

// erfi(x) = -i erf(ix)
double FADDEEVA_RE(erfi)(double x)
{
  return x*x > 720 ? (x > 0 ? Inf : -Inf)
    : exp(x*x) * FADDEEVA(w_im)(x);
}

// erfc(x) = 1 - erf(x)
double FADDEEVA_RE(erfc)(double x)
{
#if !defined(__cplusplus)
  return erfc(x); // C99 supplies erfc in math.h
#elif (__cplusplus >= 201103L) || defined(HAVE_ERFC)
  return ::erfc(x); // C++11 supplies std::erfc in cmath
#else
  if (x*x > 750) // underflow
    return (x >= 0 ? 0.0 : 2.0);
  return x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 
    : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x);
#endif
}

// erfc(z) = 1 - erf(z)
cmplx FADDEEVA(erfc)(cmplx z, double relerr)
{
  double x = creal(z), y = cimag(z);

  if (x == 0.)
    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 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
  double mIm_z2 = -2*x*y; // Im(-z^2)
  if (mRe_z2 < -750) // underflow
    return (x >= 0 ? 0.0 : 2.0);

  if (x >= 0)
    return cexp(C(mRe_z2, mIm_z2))
      * FADDEEVA(w)(C(-y,x), relerr);
  else
    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)
double FADDEEVA_RE(Dawson)(double x)
{
  const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
  return spi2 * FADDEEVA(w_im)(x);
}

// compute Dawson(z) = sqrt(pi)/2  *  exp(-z^2) * erfi(z)
cmplx FADDEEVA(Dawson)(cmplx z, double relerr)
{
  const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
  double x = creal(z), y = cimag(z);

  // handle axes separately for speed & proper handling of x or y = Inf or NaN
  if (y == 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 C(x, // preserve sign of 0
               y * (1.
                    + y2 * (0.6666666666666666666666666666666666666667
                            + y2 * 0.26666666666666666666666666666666666667)));
    }
    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 = 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
     problems from multiplying exponentially large and small quantities. */
  if (y >= 0) {
    if (y < 5e-3) {
      if (fabs(x) < 5e-3)
        goto taylor;
      else if (fabs(mIm_z2) < 5e-3)
        goto taylor_realaxis;
    }
    cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr);
    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;
      else if (fabs(mIm_z2) < 5e-3)
        goto taylor_realaxis;
    }
    else if (isnan(y))
      return C(x == 0 ? 0 : NaN, NaN);
    cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2);
    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));

  /* for small |y| and small |xy|, 
     use Taylor series to avoid cancellation inaccuracy:
       dawson(x + iy)
        = D + y^2 (D + x - 2Dx^2)
            + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
        + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
              + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
     where D = dawson(x) 

     However, for large |x|, 2Dx -> 1 which gives cancellation problems in
     this series (many of the leading terms cancel).  So, for large |x|,
     we need to substitute a continued-fraction expansion for D.

        dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))

     The 6 terms shown here seems to be the minimum needed to be
     accurate as soon as the simpler Taylor expansion above starts
     breaking down.  Using this 6-term expansion, factoring out the
     denominator, and simplifying with Maple, we obtain:

      Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x
        = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
      Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y
        = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4

     Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction
     expansion for the real part, and a 2-term expansion for the imaginary
     part.  (This avoids overflow problems for huge |x|.)  This yields:
     
     Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
     Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)

 */
 taylor_realaxis:
  {
    double x2 = x*x;
    if (x2 > 1600) { // |x| > 40
      double y2 = y*y;
      if (x2 > 25e14) {// |x| > 5e7
        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)))) *
        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 y2 = y*y;
      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)))));
    }
  }
}

/////////////////////////////////////////////////////////////////////////

// return sinc(x) = sin(x)/x, given both x and sin(x) 
// [since we only use this in cases where sin(x) has already been computed]
static inline double sinc(double x, double sinx) { 
  return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x; 
}

// 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)));
}

static inline double sqr(double x) { return x*x; }

// precomputed table of expa2n2[n-1] = exp(-a2*n*n)
// for double-precision a2 = 0.26865... in FADDEEVA(w), below.
static const double expa2n2[] = {
  7.64405281671221563e-01,
  3.41424527166548425e-01,
  8.91072646929412548e-02,
  1.35887299055460086e-02,
  1.21085455253437481e-03,
  6.30452613933449404e-05,
  1.91805156577114683e-06,
  3.40969447714832381e-08,
  3.54175089099469393e-10,
  2.14965079583260682e-12,
  7.62368911833724354e-15,
  1.57982797110681093e-17,
  1.91294189103582677e-20,
  1.35344656764205340e-23,
  5.59535712428588720e-27,
  1.35164257972401769e-30,
  1.90784582843501167e-34,
  1.57351920291442930e-38,
  7.58312432328032845e-43,
  2.13536275438697082e-47,
  3.51352063787195769e-52,
  3.37800830266396920e-57,
  1.89769439468301000e-62,
  6.22929926072668851e-68,
  1.19481172006938722e-73,
  1.33908181133005953e-79,
  8.76924303483223939e-86,
  3.35555576166254986e-92,
  7.50264110688173024e-99,
  9.80192200745410268e-106,
  7.48265412822268959e-113,
  3.33770122566809425e-120,
  8.69934598159861140e-128,
  1.32486951484088852e-135,
  1.17898144201315253e-143,
  6.13039120236180012e-152,
  1.86258785950822098e-160,
  3.30668408201432783e-169,
  3.43017280887946235e-178,
  2.07915397775808219e-187,
  7.36384545323984966e-197,
  1.52394760394085741e-206,
  1.84281935046532100e-216,
  1.30209553802992923e-226,
  5.37588903521080531e-237,
  1.29689584599763145e-247,
  1.82813078022866562e-258,
  1.50576355348684241e-269,
  7.24692320799294194e-281,
  2.03797051314726829e-292,
  3.34880215927873807e-304,
  0.0 // underflow (also prevents reads past array end, below)
};

/////////////////////////////////////////////////////////////////////////

cmplx FADDEEVA(w)(cmplx z, double relerr)
{
  if (creal(z) == 0.0)
    return C(FADDEEVA_RE(erfcx)(cimag(z)), 
             creal(z)); // give correct sign of 0 in cimag(w)
  else if (cimag(z) == 0)
    return C(exp(-sqr(creal(z))),
             FADDEEVA(w_im)(creal(z)));

  double a, a2, c;
  if (relerr <= DBL_EPSILON) {
    relerr = DBL_EPSILON;
    a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))
    c = 0.329973702884629072537; // (2/pi) * a;
    a2 = 0.268657157075235951582; // a^2
  }
  else {
    const double pi = 3.14159265358979323846264338327950288419716939937510582;
    if (relerr > 0.1) relerr = 0.1; // not sensible to compute < 1 digit
    a = pi / sqrt(-log(relerr*0.5));
    c = (2/pi)*a;
    a2 = a*a;
  }
  const double x = fabs(creal(z));
  const double y = cimag(z), ya = fabs(y);

  cmplx ret = 0.; // return value

  double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;

#define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z|

#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))) {
    
    /* Poppe & Wijers suggest using a number of terms
           nu = 3 + 1442 / (26*rho + 77)
       where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
       (They only use this expansion for rho >= 1, but rho a little less
        than 1 seems okay too.)
       Instead, I did my own fit to a slightly different function
       that avoids the hypotenuse calculation, using NLopt to minimize
       the sum of the squares of the errors in nu with the constraint
       that the estimated nu be >= minimum nu to attain machine precision.
       I also separate the regions where nu == 2 and nu == 1. */
    const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
    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 = 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 = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
      }
    }
    else { // compute nu(z) estimate and do general continued fraction
      const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit
      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) = i/sqrt(pi) / w:
        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(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
    }
    else
      return ret;
  }
#else // !USE_CONTINUED_FRACTION
  if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision
    const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
    double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
    // scale to avoid overflow
    if (x > ya) {
      double yax = ya / xs; 
      double denom = ispi / (xs + yax*ya);
      ret = C(denom*yax, denom);
    }
    else {
      double xya = xs / ya;
      double denom = ispi / (xya*xs + ya);
      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(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
    }
    else
      return ret;
  }
#endif // !USE_CONTINUED_FRACTION 

  /* Note: The test that seems to be suggested in the paper is x <
     sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2)
     underflows to zero and sum1,sum2,sum4 are zero.  However, long
     before this occurs, the sum1,sum2,sum4 contributions are
     negligible in double precision; I find that this happens for x >
     about 6, for all y.  On the other hand, I find that the case
     where we compute all of the sums is faster (at least with the
     precomputed expa2n2 table) until about x=10.  Furthermore, if we
     try to compute all of the sums for x > 20, I find that we
     sometimes run into numerical problems because underflow/overflow
     problems start to appear in the various coefficients of the sums,
     below.  Therefore, we use x < 10 here. */
  else if (x < 10) {
    double prod2ax = 1, prodm2ax = 1;
    double expx2;

    if (isnan(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
       exponential, and the x < 5e-4 special case is needed for accuracy. */

    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
        // 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;
        }
      }
      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;
        }
      }
    }
    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;
        }
      }
      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;
        }
      }
    }
    const double expx2erfcxy = // avoid spurious overflow for large negative y
      y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision
      ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x);
    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);
    }
    else {
      double xs = creal(z);
      const double sinxy = sin(xs*y);
      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 = 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 C(x,x);
    if (isnan(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. */
      ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
    }
    else
      ret = exp(-x*x); // not negligible in real part if y very small
#endif
    // (round instead of ceil as in original paper; note that x/a > 1 here)
    double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0
    double dx = a*n0 - x;
    sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y);
    sum5 = a*n0 * sum3;
    double exp1 = exp(4*a*dx), exp1dn = 1;
    int dn;
    for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms
      double np = n0 + dn, nm = n0 - dn;
      double tp = exp(-sqr(a*dn+dx));
      double tm = tp * (exp1dn *= exp1); // trick to get tm from tp
      tp /= (a2*(np*np) + y*y);
      tm /= (a2*(nm*nm) + y*y);
      sum3 += tp + tm;
      sum5 += a * (np * tp + nm * tm);
      if (a * (np * tp + nm * tm) < relerr * sum5) goto finish;
    }
    while (1) { // loop over n0+dn terms only (since n0-dn <= 0)
      double np = n0 + dn++;
      double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y);
      sum3 += tp;
      sum5 += a * np * tp;
      if (a * np * tp < relerr * sum5) goto finish;
    }
  }
 finish:
  return ret + C((0.5*c)*y*(sum2+sum3), 
                 (0.5*c)*copysign(sum5-sum4, creal(z)));
}

/////////////////////////////////////////////////////////////////////////

/* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by
   Steven G. Johnson, October 2012.

   This function combines a few different ideas.

   First, for x > 50, it uses a continued-fraction expansion (same as
   for the Faddeeva function, but with algebraic simplifications for z=i*x).

   Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations,
   but with two twists:

      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.

      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.

   For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
   with the usual checks for overflow etcetera.

   Performance-wise, it seems to be substantially faster than either
   the SLATEC DERFC function [or an erfcx function derived therefrom]
   or Cody's CALERF function (from netlib.org/specfun), while
   retaining near machine precision in accuracy.  */

/* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x).

   Uses a look-up table of 100 different Chebyshev polynomials
   for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
   with the help of Maple and a little shell script.   This allows
   the Chebyshev polynomials to be of significantly lower degree (about 1/4)
   compared to fitting the whole [0,1] interval with a single polynomial. */
static double erfcx_y100(double y100)
{
  switch ((int) y100) {
case 0: {
double t = 2*y100 - 1;
return 0.70878032454106438663e-3 + (0.71234091047026302958e-3 + (0.35779077297597742384e-5 + (0.17403143962587937815e-7 + (0.81710660047307788845e-10 + (0.36885022360434957634e-12 + 0.15917038551111111111e-14 * t) * t) * t) * t) * t) * t;
}
case 1: {
double t = 2*y100 - 3;
return 0.21479143208285144230e-2 + (0.72686402367379996033e-3 + (0.36843175430938995552e-5 + (0.18071841272149201685e-7 + (0.85496449296040325555e-10 + (0.38852037518534291510e-12 + 0.16868473576888888889e-14 * t) * t) * t) * t) * t) * t;
}
case 2: {
double t = 2*y100 - 5;
return 0.36165255935630175090e-2 + (0.74182092323555510862e-3 + (0.37948319957528242260e-5 + (0.18771627021793087350e-7 + (0.89484715122415089123e-10 + (0.40935858517772440862e-12 + 0.17872061464888888889e-14 * t) * t) * t) * t) * t) * t;
}
case 3: {
double t = 2*y100 - 7;
return 0.51154983860031979264e-2 + (0.75722840734791660540e-3 + (0.39096425726735703941e-5 + (0.19504168704300468210e-7 + (0.93687503063178993915e-10 + (0.43143925959079664747e-12 + 0.18939926435555555556e-14 * t) * t) * t) * t) * t) * t;
}
case 4: {
double t = 2*y100 - 9;
return 0.66457513172673049824e-2 + (0.77310406054447454920e-3 + (0.40289510589399439385e-5 + (0.20271233238288381092e-7 + (0.98117631321709100264e-10 + (0.45484207406017752971e-12 + 0.20076352213333333333e-14 * t) * t) * t) * t) * t) * t;
}
case 5: {
double t = 2*y100 - 11;
return 0.82082389970241207883e-2 + (0.78946629611881710721e-3 + (0.41529701552622656574e-5 + (0.21074693344544655714e-7 + (0.10278874108587317989e-9 + (0.47965201390613339638e-12 + 0.21285907413333333333e-14 * t) * t) * t) * t) * t) * t;
}
case 6: {
double t = 2*y100 - 13;
return 0.98039537275352193165e-2 + (0.80633440108342840956e-3 + (0.42819241329736982942e-5 + (0.21916534346907168612e-7 + (0.10771535136565470914e-9 + (0.50595972623692822410e-12 + 0.22573462684444444444e-14 * t) * t) * t) * t) * t) * t;
}
case 7: {
double t = 2*y100 - 15;
return 0.11433927298290302370e-1 + (0.82372858383196561209e-3 + (0.44160495311765438816e-5 + (0.22798861426211986056e-7 + (0.11291291745879239736e-9 + (0.53386189365816880454e-12 + 0.23944209546666666667e-14 * t) * t) * t) * t) * t) * t;
}
case 8: {
double t = 2*y100 - 17;
return 0.13099232878814653979e-1 + (0.84167002467906968214e-3 + (0.45555958988457506002e-5 + (0.23723907357214175198e-7 + (0.11839789326602695603e-9 + (0.56346163067550237877e-12 + 0.25403679644444444444e-14 * t) * t) * t) * t) * t) * t;
}
case 9: {
double t = 2*y100 - 19;
return 0.14800987015587535621e-1 + (0.86018092946345943214e-3 + (0.47008265848816866105e-5 + (0.24694040760197315333e-7 + (0.12418779768752299093e-9 + (0.59486890370320261949e-12 + 0.26957764568888888889e-14 * t) * t) * t) * t) * t) * t;
}
case 10: {
double t = 2*y100 - 21;
return 0.16540351739394069380e-1 + (0.87928458641241463952e-3 + (0.48520195793001753903e-5 + (0.25711774900881709176e-7 + (0.13030128534230822419e-9 + (0.62820097586874779402e-12 + 0.28612737351111111111e-14 * t) * t) * t) * t) * t) * t;
}
case 11: {
double t = 2*y100 - 23;
return 0.18318536789842392647e-1 + (0.89900542647891721692e-3 + (0.50094684089553365810e-5 + (0.26779777074218070482e-7 + (0.13675822186304615566e-9 + (0.66358287745352705725e-12 + 0.30375273884444444444e-14 * t) * t) * t) * t) * t) * t;
}
case 12: {
double t = 2*y100 - 25;
return 0.20136801964214276775e-1 + (0.91936908737673676012e-3 + (0.51734830914104276820e-5 + (0.27900878609710432673e-7 + (0.14357976402809042257e-9 + (0.70114790311043728387e-12 + 0.32252476000000000000e-14 * t) * t) * t) * t) * t) * t;
}
case 13: {
double t = 2*y100 - 27;
return 0.21996459598282740954e-1 + (0.94040248155366777784e-3 + (0.53443911508041164739e-5 + (0.29078085538049374673e-7 + (0.15078844500329731137e-9 + (0.74103813647499204269e-12 + 0.34251892320000000000e-14 * t) * t) * t) * t) * t) * t;
}
case 14: {
double t = 2*y100 - 29;
return 0.23898877187226319502e-1 + (0.96213386835900177540e-3 + (0.55225386998049012752e-5 + (0.30314589961047687059e-7 + (0.15840826497296335264e-9 + (0.78340500472414454395e-12 + 0.36381553564444444445e-14 * t) * t) * t) * t) * t) * t;
}
case 15: {
double t = 2*y100 - 31;
return 0.25845480155298518485e-1 + (0.98459293067820123389e-3 + (0.57082915920051843672e-5 + (0.31613782169164830118e-7 + (0.16646478745529630813e-9 + (0.82840985928785407942e-12 + 0.38649975768888888890e-14 * t) * t) * t) * t) * t) * t;
}
case 16: {
double t = 2*y100 - 33;
return 0.27837754783474696598e-1 + (0.10078108563256892757e-2 + (0.59020366493792212221e-5 + (0.32979263553246520417e-7 + (0.17498524159268458073e-9 + (0.87622459124842525110e-12 + 0.41066206488888888890e-14 * t) * t) * t) * t) * t) * t;
}
case 17: {
double t = 2*y100 - 35;
return 0.29877251304899307550e-1 + (0.10318204245057349310e-2 + (0.61041829697162055093e-5 + (0.34414860359542720579e-7 + (0.18399863072934089607e-9 + (0.92703227366365046533e-12 + 0.43639844053333333334e-14 * t) * t) * t) * t) * t) * t;
}
case 18: {
double t = 2*y100 - 37;
return 0.31965587178596443475e-1 + (0.10566560976716574401e-2 + (0.63151633192414586770e-5 + (0.35924638339521924242e-7 + (0.19353584758781174038e-9 + (0.98102783859889264382e-12 + 0.46381060817777777779e-14 * t) * t) * t) * t) * t) * t;
}
case 19: {
double t = 2*y100 - 39;
return 0.34104450552588334840e-1 + (0.10823541191350532574e-2 + (0.65354356159553934436e-5 + (0.37512918348533521149e-7 + (0.20362979635817883229e-9 + (0.10384187833037282363e-11 + 0.49300625262222222221e-14 * t) * t) * t) * t) * t) * t;
}
case 20: {
double t = 2*y100 - 41;
return 0.36295603928292425716e-1 + (0.11089526167995268200e-2 + (0.67654845095518363577e-5 + (0.39184292949913591646e-7 + (0.21431552202133775150e-9 + (0.10994259106646731797e-11 + 0.52409949102222222221e-14 * t) * t) * t) * t) * t) * t;
}
case 21: {
double t = 2*y100 - 43;
return 0.38540888038840509795e-1 + (0.11364917134175420009e-2 + (0.70058230641246312003e-5 + (0.40943644083718586939e-7 + (0.22563034723692881631e-9 + (0.11642841011361992885e-11 + 0.55721092871111111110e-14 * t) * t) * t) * t) * t) * t;
}
case 22: {
double t = 2*y100 - 45;
return 0.40842225954785960651e-1 + (0.11650136437945673891e-2 + (0.72569945502343006619e-5 + (0.42796161861855042273e-7 + (0.23761401711005024162e-9 + (0.12332431172381557035e-11 + 0.59246802364444444445e-14 * t) * t) * t) * t) * t) * t;
}
case 23: {
double t = 2*y100 - 47;
return 0.43201627431540222422e-1 + (0.11945628793917272199e-2 + (0.75195743532849206263e-5 + (0.44747364553960993492e-7 + (0.25030885216472953674e-9 + (0.13065684400300476484e-11 + 0.63000532853333333334e-14 * t) * t) * t) * t) * t) * t;
}
case 24: {
double t = 2*y100 - 49;
return 0.45621193513810471438e-1 + (0.12251862608067529503e-2 + (0.77941720055551920319e-5 + (0.46803119830954460212e-7 + (0.26375990983978426273e-9 + (0.13845421370977119765e-11 + 0.66996477404444444445e-14 * t) * t) * t) * t) * t) * t;
}
case 25: {
double t = 2*y100 - 51;
return 0.48103121413299865517e-1 + (0.12569331386432195113e-2 + (0.80814333496367673980e-5 + (0.48969667335682018324e-7 + (0.27801515481905748484e-9 + (0.14674637611609884208e-11 + 0.71249589351111111110e-14 * t) * t) * t) * t) * t) * t;
}
case 26: {
double t = 2*y100 - 53;
return 0.50649709676983338501e-1 + (0.12898555233099055810e-2 + (0.83820428414568799654e-5 + (0.51253642652551838659e-7 + (0.29312563849675507232e-9 + (0.15556512782814827846e-11 + 0.75775607822222222221e-14 * t) * t) * t) * t) * t) * t;
}
case 27: {
double t = 2*y100 - 55;
return 0.53263363664388864181e-1 + (0.13240082443256975769e-2 + (0.86967260015007658418e-5 + (0.53662102750396795566e-7 + (0.30914568786634796807e-9 + (0.16494420240828493176e-11 + 0.80591079644444444445e-14 * t) * t) * t) * t) * t) * t;
}
case 28: {
double t = 2*y100 - 57;
return 0.55946601353500013794e-1 + (0.13594491197408190706e-2 + (0.90262520233016380987e-5 + (0.56202552975056695376e-7 + (0.32613310410503135996e-9 + (0.17491936862246367398e-11 + 0.85713381688888888890e-14 * t) * t) * t) * t) * t) * t;
}
case 29: {
double t = 2*y100 - 59;
return 0.58702059496154081813e-1 + (0.13962391363223647892e-2 + (0.93714365487312784270e-5 + (0.58882975670265286526e-7 + (0.34414937110591753387e-9 + (0.18552853109751857859e-11 + 0.91160736711111111110e-14 * t) * t) * t) * t) * t) * t;
}
case 30: {
double t = 2*y100 - 61;
return 0.61532500145144778048e-1 + (0.14344426411912015247e-2 + (0.97331446201016809696e-5 + (0.61711860507347175097e-7 + (0.36325987418295300221e-9 + (0.19681183310134518232e-11 + 0.96952238400000000000e-14 * t) * t) * t) * t) * t) * t;
}
case 31: {
double t = 2*y100 - 63;
return 0.64440817576653297993e-1 + (0.14741275456383131151e-2 + (0.10112293819576437838e-4 + (0.64698236605933246196e-7 + (0.38353412915303665586e-9 + (0.20881176114385120186e-11 + 0.10310784480000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 32: {
double t = 2*y100 - 65;
return 0.67430045633130393282e-1 + (0.15153655418916540370e-2 + (0.10509857606888328667e-4 + (0.67851706529363332855e-7 + (0.40504602194811140006e-9 + (0.22157325110542534469e-11 + 0.10964842115555555556e-13 * t) * t) * t) * t) * t) * t;
}
case 33: {
double t = 2*y100 - 67;
return 0.70503365513338850709e-1 + (0.15582323336495709827e-2 + (0.10926868866865231089e-4 + (0.71182482239613507542e-7 + (0.42787405890153386710e-9 + (0.23514379522274416437e-11 + 0.11659571751111111111e-13 * t) * t) * t) * t) * t) * t;
}
case 34: {
double t = 2*y100 - 69;
return 0.73664114037944596353e-1 + (0.16028078812438820413e-2 + (0.11364423678778207991e-4 + (0.74701423097423182009e-7 + (0.45210162777476488324e-9 + (0.24957355004088569134e-11 + 0.12397238257777777778e-13 * t) * t) * t) * t) * t) * t;
}
case 35: {
double t = 2*y100 - 71;
return 0.76915792420819562379e-1 + (0.16491766623447889354e-2 + (0.11823685320041302169e-4 + (0.78420075993781544386e-7 + (0.47781726956916478925e-9 + (0.26491544403815724749e-11 + 0.13180196462222222222e-13 * t) * t) * t) * t) * t) * t;
}
case 36: {
double t = 2*y100 - 73;
return 0.80262075578094612819e-1 + (0.16974279491709504117e-2 + (0.12305888517309891674e-4 + (0.82350717698979042290e-7 + (0.50511496109857113929e-9 + (0.28122528497626897696e-11 + 0.14010889635555555556e-13 * t) * t) * t) * t) * t) * t;
}
case 37: {
double t = 2*y100 - 75;
return 0.83706822008980357446e-1 + (0.17476561032212656962e-2 + (0.12812343958540763368e-4 + (0.86506399515036435592e-7 + (0.53409440823869467453e-9 + (0.29856186620887555043e-11 + 0.14891851591111111111e-13 * t) * t) * t) * t) * t) * t;
}
case 38: {
double t = 2*y100 - 77;
return 0.87254084284461718231e-1 + (0.17999608886001962327e-2 + (0.13344443080089492218e-4 + (0.90900994316429008631e-7 + (0.56486134972616465316e-9 + (0.31698707080033956934e-11 + 0.15825697795555555556e-13 * t) * t) * t) * t) * t) * t;
}
case 39: {
double t = 2*y100 - 79;
return 0.90908120182172748487e-1 + (0.18544478050657699758e-2 + (0.13903663143426120077e-4 + (0.95549246062549906177e-7 + (0.59752787125242054315e-9 + (0.33656597366099099413e-11 + 0.16815130613333333333e-13 * t) * t) * t) * t) * t) * t;
}
case 40: {
double t = 2*y100 - 81;
return 0.94673404508075481121e-1 + (0.19112284419887303347e-2 + (0.14491572616545004930e-4 + (0.10046682186333613697e-6 + (0.63221272959791000515e-9 + (0.35736693975589130818e-11 + 0.17862931591111111111e-13 * t) * t) * t) * t) * t) * t;
}
case 41: {
double t = 2*y100 - 83;
return 0.98554641648004456555e-1 + (0.19704208544725622126e-2 + (0.15109836875625443935e-4 + (0.10567036667675984067e-6 + (0.66904168640019354565e-9 + (0.37946171850824333014e-11 + 0.18971959040000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 42: {
double t = 2*y100 - 85;
return 0.10255677889470089531e0 + (0.20321499629472857418e-2 + (0.15760224242962179564e-4 + (0.11117756071353507391e-6 + (0.70814785110097658502e-9 + (0.40292553276632563925e-11 + 0.20145143075555555556e-13 * t) * t) * t) * t) * t) * t;
}
case 43: {
double t = 2*y100 - 87;
return 0.10668502059865093318e0 + (0.20965479776148731610e-2 + (0.16444612377624983565e-4 + (0.11700717962026152749e-6 + (0.74967203250938418991e-9 + (0.42783716186085922176e-11 + 0.21385479360000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 44: {
double t = 2*y100 - 89;
return 0.11094484319386444474e0 + (0.21637548491908170841e-2 + (0.17164995035719657111e-4 + (0.12317915750735938089e-6 + (0.79376309831499633734e-9 + (0.45427901763106353914e-11 + 0.22696025653333333333e-13 * t) * t) * t) * t) * t) * t;
}
case 45: {
double t = 2*y100 - 91;
return 0.11534201115268804714e0 + (0.22339187474546420375e-2 + (0.17923489217504226813e-4 + (0.12971465288245997681e-6 + (0.84057834180389073587e-9 + (0.48233721206418027227e-11 + 0.24079890062222222222e-13 * t) * t) * t) * t) * t) * t;
}
case 46: {
double t = 2*y100 - 93;
return 0.11988259392684094740e0 + (0.23071965691918689601e-2 + (0.18722342718958935446e-4 + (0.13663611754337957520e-6 + (0.89028385488493287005e-9 + (0.51210161569225846701e-11 + 0.25540227111111111111e-13 * t) * t) * t) * t) * t) * t;
}
case 47: {
double t = 2*y100 - 95;
return 0.12457298393509812907e0 + (0.23837544771809575380e-2 + (0.19563942105711612475e-4 + (0.14396736847739470782e-6 + (0.94305490646459247016e-9 + (0.54366590583134218096e-11 + 0.27080225920000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 48: {
double t = 2*y100 - 97;
return 0.12941991566142438816e0 + (0.24637684719508859484e-2 + (0.20450821127475879816e-4 + (0.15173366280523906622e-6 + (0.99907632506389027739e-9 + (0.57712760311351625221e-11 + 0.28703099555555555556e-13 * t) * t) * t) * t) * t) * t;
}
case 49: {
double t = 2*y100 - 99;
return 0.13443048593088696613e0 + (0.25474249981080823877e-2 + (0.21385669591362915223e-4 + (0.15996177579900443030e-6 + (0.10585428844575134013e-8 + (0.61258809536787882989e-11 + 0.30412080142222222222e-13 * t) * t) * t) * t) * t) * t;
}
case 50: {
double t = 2*y100 - 101;
return 0.13961217543434561353e0 + (0.26349215871051761416e-2 + (0.22371342712572567744e-4 + (0.16868008199296822247e-6 + (0.11216596910444996246e-8 + (0.65015264753090890662e-11 + 0.32210394506666666666e-13 * t) * t) * t) * t) * t) * t;
}
case 51: {
double t = 2*y100 - 103;
return 0.14497287157673800690e0 + (0.27264675383982439814e-2 + (0.23410870961050950197e-4 + (0.17791863939526376477e-6 + (0.11886425714330958106e-8 + (0.68993039665054288034e-11 + 0.34101266222222222221e-13 * t) * t) * t) * t) * t) * t;
}
case 52: {
double t = 2*y100 - 105;
return 0.15052089272774618151e0 + (0.28222846410136238008e-2 + (0.24507470422713397006e-4 + (0.18770927679626136909e-6 + (0.12597184587583370712e-8 + (0.73203433049229821618e-11 + 0.36087889048888888890e-13 * t) * t) * t) * t) * t) * t;
}
case 53: {
double t = 2*y100 - 107;
return 0.15626501395774612325e0 + (0.29226079376196624949e-2 + (0.25664553693768450545e-4 + (0.19808568415654461964e-6 + (0.13351257759815557897e-8 + (0.77658124891046760667e-11 + 0.38173420035555555555e-13 * t) * t) * t) * t) * t) * t;
}
case 54: {
double t = 2*y100 - 109;
return 0.16221449434620737567e0 + (0.30276865332726475672e-2 + (0.26885741326534564336e-4 + (0.20908350604346384143e-6 + (0.14151148144240728728e-8 + (0.82369170665974313027e-11 + 0.40360957457777777779e-13 * t) * t) * t) * t) * t) * t;
}
case 55: {
double t = 2*y100 - 111;
return 0.16837910595412130659e0 + (0.31377844510793082301e-2 + (0.28174873844911175026e-4 + (0.22074043807045782387e-6 + (0.14999481055996090039e-8 + (0.87348993661930809254e-11 + 0.42653528977777777779e-13 * t) * t) * t) * t) * t) * t;
}
case 56: {
double t = 2*y100 - 113;
return 0.17476916455659369953e0 + (0.32531815370903068316e-2 + (0.29536024347344364074e-4 + (0.23309632627767074202e-6 + (0.15899007843582444846e-8 + (0.92610375235427359475e-11 + 0.45054073102222222221e-13 * t) * t) * t) * t) * t) * t;
}
case 57: {
double t = 2*y100 - 115;
return 0.18139556223643701364e0 + (0.33741744168096996041e-2 + (0.30973511714709500836e-4 + (0.24619326937592290996e-6 + (0.16852609412267750744e-8 + (0.98166442942854895573e-11 + 0.47565418097777777779e-13 * t) * t) * t) * t) * t) * t;
}
case 58: {
double t = 2*y100 - 117;
return 0.18826980194443664549e0 + (0.35010775057740317997e-2 + (0.32491914440014267480e-4 + (0.26007572375886319028e-6 + (0.17863299617388376116e-8 + (0.10403065638343878679e-10 + 0.50190265831111111110e-13 * t) * t) * t) * t) * t) * t;
}
case 59: {
double t = 2*y100 - 119;
return 0.19540403413693967350e0 + (0.36342240767211326315e-2 + (0.34096085096200907289e-4 + (0.27479061117017637474e-6 + (0.18934228504790032826e-8 + (0.11021679075323598664e-10 + 0.52931171733333333334e-13 * t) * t) * t) * t) * t) * t;
}
case 60: {
double t = 2*y100 - 121;
return 0.20281109560651886959e0 + (0.37739673859323597060e-2 + (0.35791165457592409054e-4 + (0.29038742889416172404e-6 + (0.20068685374849001770e-8 + (0.11673891799578381999e-10 + 0.55790523093333333334e-13 * t) * t) * t) * t) * t) * t;
}
case 61: {
double t = 2*y100 - 123;
return 0.21050455062669334978e0 + (0.39206818613925652425e-2 + (0.37582602289680101704e-4 + (0.30691836231886877385e-6 + (0.21270101645763677824e-8 + (0.12361138551062899455e-10 + 0.58770520160000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 62: {
double t = 2*y100 - 125;
return 0.21849873453703332479e0 + (0.40747643554689586041e-2 + (0.39476163820986711501e-4 + (0.32443839970139918836e-6 + (0.22542053491518680200e-8 + (0.13084879235290858490e-10 + 0.61873153262222222221e-13 * t) * t) * t) * t) * t) * t;
}
case 63: {
double t = 2*y100 - 127;
return 0.22680879990043229327e0 + (0.42366354648628516935e-2 + (0.41477956909656896779e-4 + (0.34300544894502810002e-6 + (0.23888264229264067658e-8 + (0.13846596292818514601e-10 + 0.65100183751111111110e-13 * t) * t) * t) * t) * t) * t;
}
case 64: {
double t = 2*y100 - 129;
return 0.23545076536988703937e0 + (0.44067409206365170888e-2 + (0.43594444916224700881e-4 + (0.36268045617760415178e-6 + (0.25312606430853202748e-8 + (0.14647791812837903061e-10 + 0.68453122631111111110e-13 * t) * t) * t) * t) * t) * t;
}
case 65: {
double t = 2*y100 - 131;
return 0.24444156740777432838e0 + (0.45855530511605787178e-2 + (0.45832466292683085475e-4 + (0.38352752590033030472e-6 + (0.26819103733055603460e-8 + (0.15489984390884756993e-10 + 0.71933206364444444445e-13 * t) * t) * t) * t) * t) * t;
}
case 66: {
double t = 2*y100 - 133;
return 0.25379911500634264643e0 + (0.47735723208650032167e-2 + (0.48199253896534185372e-4 + (0.40561404245564732314e-6 + (0.28411932320871165585e-8 + (0.16374705736458320149e-10 + 0.75541379822222222221e-13 * t) * t) * t) * t) * t) * t;
}
case 67: {
double t = 2*y100 - 135;
return 0.26354234756393613032e0 + (0.49713289477083781266e-2 + (0.50702455036930367504e-4 + (0.42901079254268185722e-6 + (0.30095422058900481753e-8 + (0.17303497025347342498e-10 + 0.79278273368888888890e-13 * t) * t) * t) * t) * t) * t;
}
case 68: {
double t = 2*y100 - 137;
return 0.27369129607732343398e0 + (0.51793846023052643767e-2 + (0.53350152258326602629e-4 + (0.45379208848865015485e-6 + (0.31874057245814381257e-8 + (0.18277905010245111046e-10 + 0.83144182364444444445e-13 * t) * t) * t) * t) * t) * t;
}
case 69: {
double t = 2*y100 - 139;
return 0.28426714781640316172e0 + (0.53983341916695141966e-2 + (0.56150884865255810638e-4 + (0.48003589196494734238e-6 + (0.33752476967570796349e-8 + (0.19299477888083469086e-10 + 0.87139049137777777779e-13 * t) * t) * t) * t) * t) * t;
}
case 70: {
double t = 2*y100 - 141;
return 0.29529231465348519920e0 + (0.56288077305420795663e-2 + (0.59113671189913307427e-4 + (0.50782393781744840482e-6 + (0.35735475025851713168e-8 + (0.20369760937017070382e-10 + 0.91262442613333333334e-13 * t) * t) * t) * t) * t) * t;
}
case 71: {
double t = 2*y100 - 143;
return 0.30679050522528838613e0 + (0.58714723032745403331e-2 + (0.62248031602197686791e-4 + (0.53724185766200945789e-6 + (0.37827999418960232678e-8 + (0.21490291930444538307e-10 + 0.95513539182222222221e-13 * t) * t) * t) * t) * t) * t;
}
case 72: {
double t = 2*y100 - 145;
return 0.31878680111173319425e0 + (0.61270341192339103514e-2 + (0.65564012259707640976e-4 + (0.56837930287837738996e-6 + (0.40035151353392378882e-8 + (0.22662596341239294792e-10 + 0.99891109760000000000e-13 * t) * t) * t) * t) * t) * t;
}
case 73: {
double t = 2*y100 - 147;
return 0.33130773722152622027e0 + (0.63962406646798080903e-2 + (0.69072209592942396666e-4 + (0.60133006661885941812e-6 + (0.42362183765883466691e-8 + (0.23888182347073698382e-10 + 0.10439349811555555556e-12 * t) * t) * t) * t) * t) * t;
}
case 74: {
double t = 2*y100 - 149;
return 0.34438138658041336523e0 + (0.66798829540414007258e-2 + (0.72783795518603561144e-4 + (0.63619220443228800680e-6 + (0.44814499336514453364e-8 + (0.25168535651285475274e-10 + 0.10901861383111111111e-12 * t) * t) * t) * t) * t) * t;
}
case 75: {
double t = 2*y100 - 151;
return 0.35803744972380175583e0 + (0.69787978834882685031e-2 + (0.76710543371454822497e-4 + (0.67306815308917386747e-6 + (0.47397647975845228205e-8 + (0.26505114141143050509e-10 + 0.11376390933333333333e-12 * t) * t) * t) * t) * t) * t;
}
case 76: {
double t = 2*y100 - 153;
return 0.37230734890119724188e0 + (0.72938706896461381003e-2 + (0.80864854542670714092e-4 + (0.71206484718062688779e-6 + (0.50117323769745883805e-8 + (0.27899342394100074165e-10 + 0.11862637614222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 77: {
double t = 2*y100 - 155;
return 0.38722432730555448223e0 + (0.76260375162549802745e-2 + (0.85259785810004603848e-4 + (0.75329383305171327677e-6 + (0.52979361368388119355e-8 + (0.29352606054164086709e-10 + 0.12360253370666666667e-12 * t) * t) * t) * t) * t) * t;
}
case 78: {
double t = 2*y100 - 157;
return 0.40282355354616940667e0 + (0.79762880915029728079e-2 + (0.89909077342438246452e-4 + (0.79687137961956194579e-6 + (0.55989731807360403195e-8 + (0.30866246101464869050e-10 + 0.12868841946666666667e-12 * t) * t) * t) * t) * t) * t;
}
case 79: {
double t = 2*y100 - 159;
return 0.41914223158913787649e0 + (0.83456685186950463538e-2 + (0.94827181359250161335e-4 + (0.84291858561783141014e-6 + (0.59154537751083485684e-8 + (0.32441553034347469291e-10 + 0.13387957943111111111e-12 * t) * t) * t) * t) * t) * t;
}
case 80: {
double t = 2*y100 - 161;
return 0.43621971639463786896e0 + (0.87352841828289495773e-2 + (0.10002929142066799966e-3 + (0.89156148280219880024e-6 + (0.62480008150788597147e-8 + (0.34079760983458878910e-10 + 0.13917107176888888889e-12 * t) * t) * t) * t) * t) * t;
}
case 81: {
double t = 2*y100 - 163;
return 0.45409763548534330981e0 + (0.91463027755548240654e-2 + (0.10553137232446167258e-3 + (0.94293113464638623798e-6 + (0.65972492312219959885e-8 + (0.35782041795476563662e-10 + 0.14455745872000000000e-12 * t) * t) * t) * t) * t) * t;
}
case 82: {
double t = 2*y100 - 165;
return 0.47282001668512331468e0 + (0.95799574408860463394e-2 + (0.11135019058000067469e-3 + (0.99716373005509038080e-6 + (0.69638453369956970347e-8 + (0.37549499088161345850e-10 + 0.15003280712888888889e-12 * t) * t) * t) * t) * t) * t;
}
case 83: {
double t = 2*y100 - 167;
return 0.49243342227179841649e0 + (0.10037550043909497071e-1 + (0.11750334542845234952e-3 + (0.10544006716188967172e-5 + (0.73484461168242224872e-8 + (0.39383162326435752965e-10 + 0.15559069118222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 84: {
double t = 2*y100 - 169;
return 0.51298708979209258326e0 + (0.10520454564612427224e-1 + (0.12400930037494996655e-3 + (0.11147886579371265246e-5 + (0.77517184550568711454e-8 + (0.41283980931872622611e-10 + 0.16122419680000000000e-12 * t) * t) * t) * t) * t) * t;
}
case 85: {
double t = 2*y100 - 171;
return 0.53453307979101369843e0 + (0.11030120618800726938e-1 + (0.13088741519572269581e-3 + (0.11784797595374515432e-5 + (0.81743383063044825400e-8 + (0.43252818449517081051e-10 + 0.16692592640000000000e-12 * t) * t) * t) * t) * t) * t;
}
case 86: {
double t = 2*y100 - 173;
return 0.55712643071169299478e0 + (0.11568077107929735233e-1 + (0.13815797838036651289e-3 + (0.12456314879260904558e-5 + (0.86169898078969313597e-8 + (0.45290446811539652525e-10 + 0.17268801084444444444e-12 * t) * t) * t) * t) * t) * t;
}
case 87: {
double t = 2*y100 - 175;
return 0.58082532122519320968e0 + (0.12135935999503877077e-1 + (0.14584223996665838559e-3 + (0.13164068573095710742e-5 + (0.90803643355106020163e-8 + (0.47397540713124619155e-10 + 0.17850211608888888889e-12 * t) * t) * t) * t) * t) * t;
}
case 88: {
double t = 2*y100 - 177;
return 0.60569124025293375554e0 + (0.12735396239525550361e-1 + (0.15396244472258863344e-3 + (0.13909744385382818253e-5 + (0.95651595032306228245e-8 + (0.49574672127669041550e-10 + 0.18435945564444444444e-12 * t) * t) * t) * t) * t) * t;
}
case 89: {
double t = 2*y100 - 179;
return 0.63178916494715716894e0 + (0.13368247798287030927e-1 + (0.16254186562762076141e-3 + (0.14695084048334056083e-5 + (0.10072078109604152350e-7 + (0.51822304995680707483e-10 + 0.19025081422222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 90: {
double t = 2*y100 - 181;
return 0.65918774689725319200e0 + (0.14036375850601992063e-1 + (0.17160483760259706354e-3 + (0.15521885688723188371e-5 + (0.10601827031535280590e-7 + (0.54140790105837520499e-10 + 0.19616655146666666667e-12 * t) * t) * t) * t) * t) * t;
}
case 91: {
double t = 2*y100 - 183;
return 0.68795950683174433822e0 + (0.14741765091365869084e-1 + (0.18117679143520433835e-3 + (0.16392004108230585213e-5 + (0.11155116068018043001e-7 + (0.56530360194925690374e-10 + 0.20209663662222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 92: {
double t = 2*y100 - 185;
return 0.71818103808729967036e0 + (0.15486504187117112279e-1 + (0.19128428784550923217e-3 + (0.17307350969359975848e-5 + (0.11732656736113607751e-7 + (0.58991125287563833603e-10 + 0.20803065333333333333e-12 * t) * t) * t) * t) * t) * t;
}
case 93: {
double t = 2*y100 - 187;
return 0.74993321911726254661e0 + (0.16272790364044783382e-1 + (0.20195505163377912645e-3 + (0.18269894883203346953e-5 + (0.12335161021630225535e-7 + (0.61523068312169087227e-10 + 0.21395783431111111111e-12 * t) * t) * t) * t) * t) * t;
}
case 94: {
double t = 2*y100 - 189;
return 0.78330143531283492729e0 + (0.17102934132652429240e-1 + (0.21321800585063327041e-3 + (0.19281661395543913713e-5 + (0.12963340087354341574e-7 + (0.64126040998066348872e-10 + 0.21986708942222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 95: {
double t = 2*y100 - 191;
return 0.81837581041023811832e0 + (0.17979364149044223802e-1 + (0.22510330592753129006e-3 + (0.20344732868018175389e-5 + (0.13617902941839949718e-7 + (0.66799760083972474642e-10 + 0.22574701262222222222e-12 * t) * t) * t) * t) * t) * t;
}
case 96: {
double t = 2*y100 - 193;
return 0.85525144775685126237e0 + (0.18904632212547561026e-1 + (0.23764237370371255638e-3 + (0.21461248251306387979e-5 + (0.14299555071870523786e-7 + (0.69543803864694171934e-10 + 0.23158593688888888889e-12 * t) * t) * t) * t) * t) * t;
}
case 97: {
double t = 2*y100 - 195;
return 0.89402868170849933734e0 + (0.19881418399127202569e-1 + (0.25086793128395995798e-3 + (0.22633402747585233180e-5 + (0.15008997042116532283e-7 + (0.72357609075043941261e-10 + 0.23737194737777777778e-12 * t) * t) * t) * t) * t) * t;
}
case 98: {
double t = 2*y100 - 197;
return 0.93481333942870796363e0 + (0.20912536329780368893e-1 + (0.26481403465998477969e-3 + (0.23863447359754921676e-5 + (0.15746923065472184451e-7 + (0.75240468141720143653e-10 + 0.24309291271111111111e-12 * t) * t) * t) * t) * t) * t;
}
case 99: {
double t = 2*y100 - 199;
return 0.97771701335885035464e0 + (0.22000938572830479551e-1 + (0.27951610702682383001e-3 + (0.25153688325245314530e-5 + (0.16514019547822821453e-7 + (0.78191526829368231251e-10 + 0.24873652355555555556e-12 * t) * t) * t) * t) * t) * t;
}
  }
  // we only get here if y = 1, i.e. |x| < 4*eps, in which case
  // erfcx is within 1e-15 of 1..
  return 1.0;
}

double FADDEEVA_RE(erfcx)(double x)
{
  if (x >= 0) {
    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;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
                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)));
}

/////////////////////////////////////////////////////////////////////////
/* Compute a scaled Dawson integral 
            FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
   equivalent to the imaginary part w(x) for real x.

   Uses methods similar to the erfcx calculation above: continued fractions
   for large |x|, a lookup table of Chebyshev polynomials for smaller |x|,
   and finally a Taylor expansion for |x|<0.01.
   
   Steven G. Johnson, October 2012. */

/* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).

   Uses a look-up table of 100 different Chebyshev polynomials
   for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
   with the help of Maple and a little shell script.   This allows
   the Chebyshev polynomials to be of significantly lower degree (about 1/30)
   compared to fitting the whole [0,1] interval with a single polynomial. */
static double w_im_y100(double y100, double x) {
  switch ((int) y100) {
    case 0: {
      double t = 2*y100 - 1;
      return 0.28351593328822191546e-2 + (0.28494783221378400759e-2 + (0.14427470563276734183e-4 + (0.10939723080231588129e-6 + (0.92474307943275042045e-9 + (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * t) * t) * t) * t) * t;
    }
    case 1: {
      double t = 2*y100 - 3;
      return 0.85927161243940350562e-2 + (0.29085312941641339862e-2 + (0.15106783707725582090e-4 + (0.11716709978531327367e-6 + (0.10197387816021040024e-8 + (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 2: {
      double t = 2*y100 - 5;
      return 0.14471159831187703054e-1 + (0.29703978970263836210e-2 + (0.15835096760173030976e-4 + (0.12574803383199211596e-6 + (0.11278672159518415848e-8 + (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 3: {
      double t = 2*y100 - 7;
      return 0.20476320420324610618e-1 + (0.30352843012898665856e-2 + (0.16617609387003727409e-4 + (0.13525429711163116103e-6 + (0.12515095552507169013e-8 + (0.13235687543603382345e-10 + 0.15326595042666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 4: {
      double t = 2*y100 - 9;
      return 0.26614461952489004566e-1 + (0.31034189276234947088e-2 + (0.17460268109986214274e-4 + (0.14582130824485709573e-6 + (0.13935959083809746345e-8 + (0.15249438072998932900e-10 + 0.18344741882133333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 5: {
      double t = 2*y100 - 11;
      return 0.32892330248093586215e-1 + (0.31750557067975068584e-2 + (0.18369907582308672632e-4 + (0.15761063702089457882e-6 + (0.15577638230480894382e-8 + (0.17663868462699097951e-10 + (0.22126732680711111111e-12 + 0.30273474177737853668e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 6: {
      double t = 2*y100 - 13;
      return 0.39317207681134336024e-1 + (0.32504779701937539333e-2 + (0.19354426046513400534e-4 + (0.17081646971321290539e-6 + (0.17485733959327106250e-8 + (0.20593687304921961410e-10 + (0.26917401949155555556e-12 + 0.38562123837725712270e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 7: {
      double t = 2*y100 - 15;
      return 0.45896976511367738235e-1 + (0.33300031273110976165e-2 + (0.20423005398039037313e-4 + (0.18567412470376467303e-6 + (0.19718038363586588213e-8 + (0.24175006536781219807e-10 + (0.33059982791466666666e-12 + 0.49756574284439426165e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 8: {
      double t = 2*y100 - 17;
      return 0.52640192524848962855e-1 + (0.34139883358846720806e-2 + (0.21586390240603337337e-4 + (0.20247136501568904646e-6 + (0.22348696948197102935e-8 + (0.28597516301950162548e-10 + (0.41045502119111111110e-12 + 0.65151614515238361946e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 9: {
      double t = 2*y100 - 19;
      return 0.59556171228656770456e-1 + (0.35028374386648914444e-2 + (0.22857246150998562824e-4 + (0.22156372146525190679e-6 + (0.25474171590893813583e-8 + (0.34122390890697400584e-10 + (0.51593189879111111110e-12 + 0.86775076853908006938e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 10: {
      double t = 2*y100 - 21;
      return 0.66655089485108212551e-1 + (0.35970095381271285568e-2 + (0.24250626164318672928e-4 + (0.24339561521785040536e-6 + (0.29221990406518411415e-8 + (0.41117013527967776467e-10 + (0.65786450716444444445e-12 + 0.11791885745450623331e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 11: {
      double t = 2*y100 - 23;
      return 0.73948106345519174661e-1 + (0.36970297216569341748e-2 + (0.25784588137312868792e-4 + (0.26853012002366752770e-6 + (0.33763958861206729592e-8 + (0.50111549981376976397e-10 + (0.85313857496888888890e-12 + 0.16417079927706899860e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 12: {
      double t = 2*y100 - 25;
      return 0.81447508065002963203e-1 + (0.38035026606492705117e-2 + (0.27481027572231851896e-4 + (0.29769200731832331364e-6 + (0.39336816287457655076e-8 + (0.61895471132038157624e-10 + (0.11292303213511111111e-11 + 0.23558532213703884304e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 13: {
      double t = 2*y100 - 27;
      return 0.89166884027582716628e-1 + (0.39171301322438946014e-2 + (0.29366827260422311668e-4 + (0.33183204390350724895e-6 + (0.46276006281647330524e-8 + (0.77692631378169813324e-10 + (0.15335153258844444444e-11 + 0.35183103415916026911e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 14: {
      double t = 2*y100 - 29;
      return 0.97121342888032322019e-1 + (0.40387340353207909514e-2 + (0.31475490395950776930e-4 + (0.37222714227125135042e-6 + (0.55074373178613809996e-8 + (0.99509175283990337944e-10 + (0.21552645758222222222e-11 + 0.55728651431872687605e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 15: {
      double t = 2*y100 - 31;
      return 0.10532778218603311137e0 + (0.41692873614065380607e-2 + (0.33849549774889456984e-4 + (0.42064596193692630143e-6 + (0.66494579697622432987e-8 + (0.13094103581931802337e-9 + (0.31896187409777777778e-11 + 0.97271974184476560742e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 16: {
      double t = 2*y100 - 33;
      return 0.11380523107427108222e0 + (0.43099572287871821013e-2 + (0.36544324341565929930e-4 + (0.47965044028581857764e-6 + (0.81819034238463698796e-8 + (0.17934133239549647357e-9 + (0.50956666166186293627e-11 + (0.18850487318190638010e-12 + 0.79697813173519853340e-14 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 17: {
      double t = 2*y100 - 35;
      return 0.12257529703447467345e0 + (0.44621675710026986366e-2 + (0.39634304721292440285e-4 + (0.55321553769873381819e-6 + (0.10343619428848520870e-7 + (0.26033830170470368088e-9 + (0.87743837749108025357e-11 + (0.34427092430230063401e-12 + 0.10205506615709843189e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 18: {
      double t = 2*y100 - 37;
      return 0.13166276955656699478e0 + (0.46276970481783001803e-2 + (0.43225026380496399310e-4 + (0.64799164020016902656e-6 + (0.13580082794704641782e-7 + (0.39839800853954313927e-9 + (0.14431142411840000000e-10 + 0.42193457308830027541e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 19: {
      double t = 2*y100 - 39;
      return 0.14109647869803356475e0 + (0.48088424418545347758e-2 + (0.47474504753352150205e-4 + (0.77509866468724360352e-6 + (0.18536851570794291724e-7 + (0.60146623257887570439e-9 + (0.18533978397305276318e-10 + (0.41033845938901048380e-13 - 0.46160680279304825485e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 20: {
      double t = 2*y100 - 41;
      return 0.15091057940548936603e0 + (0.50086864672004685703e-2 + (0.52622482832192230762e-4 + (0.95034664722040355212e-6 + (0.25614261331144718769e-7 + (0.80183196716888606252e-9 + (0.12282524750534352272e-10 + (-0.10531774117332273617e-11 - 0.86157181395039646412e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 21: {
      double t = 2*y100 - 43;
      return 0.16114648116017010770e0 + (0.52314661581655369795e-2 + (0.59005534545908331315e-4 + (0.11885518333915387760e-5 + (0.33975801443239949256e-7 + (0.82111547144080388610e-9 + (-0.12357674017312854138e-10 + (-0.24355112256914479176e-11 - 0.75155506863572930844e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 22: {
      double t = 2*y100 - 45;
      return 0.17185551279680451144e0 + (0.54829002967599420860e-2 + (0.67013226658738082118e-4 + (0.14897400671425088807e-5 + (0.40690283917126153701e-7 + (0.44060872913473778318e-9 + (-0.52641873433280000000e-10 - 0.30940587864543343124e-11 * t) * t) * t) * t) * t) * t) * t;
    }
    case 23: {
      double t = 2*y100 - 47;
      return 0.18310194559815257381e0 + (0.57701559375966953174e-2 + (0.76948789401735193483e-4 + (0.18227569842290822512e-5 + (0.41092208344387212276e-7 + (-0.44009499965694442143e-9 + (-0.92195414685628803451e-10 + (-0.22657389705721753299e-11 + 0.10004784908106839254e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 24: {
      double t = 2*y100 - 49;
      return 0.19496527191546630345e0 + (0.61010853144364724856e-2 + (0.88812881056342004864e-4 + (0.21180686746360261031e-5 + (0.30652145555130049203e-7 + (-0.16841328574105890409e-8 + (-0.11008129460612823934e-9 + (-0.12180794204544515779e-12 + 0.15703325634590334097e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 25: {
      double t = 2*y100 - 51;
      return 0.20754006813966575720e0 + (0.64825787724922073908e-2 + (0.10209599627522311893e-3 + (0.22785233392557600468e-5 + (0.73495224449907568402e-8 + (-0.29442705974150112783e-8 + (-0.94082603434315016546e-10 + (0.23609990400179321267e-11 + 0.14141908654269023788e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 26: {
      double t = 2*y100 - 53;
      return 0.22093185554845172146e0 + (0.69182878150187964499e-2 + (0.11568723331156335712e-3 + (0.22060577946323627739e-5 + (-0.26929730679360840096e-7 + (-0.38176506152362058013e-8 + (-0.47399503861054459243e-10 + (0.40953700187172127264e-11 + 0.69157730376118511127e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 27: {
      double t = 2*y100 - 55;
      return 0.23524827304057813918e0 + (0.74063350762008734520e-2 + (0.12796333874615790348e-3 + (0.18327267316171054273e-5 + (-0.66742910737957100098e-7 + (-0.40204740975496797870e-8 + (0.14515984139495745330e-10 + (0.44921608954536047975e-11 - 0.18583341338983776219e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 28: {
      double t = 2*y100 - 57;
      return 0.25058626331812744775e0 + (0.79377285151602061328e-2 + (0.13704268650417478346e-3 + (0.11427511739544695861e-5 + (-0.10485442447768377485e-6 + (-0.34850364756499369763e-8 + (0.72656453829502179208e-10 + (0.36195460197779299406e-11 - 0.84882136022200714710e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 29: {
      double t = 2*y100 - 59;
      return 0.26701724900280689785e0 + (0.84959936119625864274e-2 + (0.14112359443938883232e-3 + (0.17800427288596909634e-6 + (-0.13443492107643109071e-6 + (-0.23512456315677680293e-8 + (0.11245846264695936769e-9 + (0.19850501334649565404e-11 - 0.11284666134635050832e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 30: {
      double t = 2*y100 - 61;
      return 0.28457293586253654144e0 + (0.90581563892650431899e-2 + (0.13880520331140646738e-3 + (-0.97262302362522896157e-6 + (-0.15077100040254187366e-6 + (-0.88574317464577116689e-9 + (0.12760311125637474581e-9 + (0.20155151018282695055e-12 - 0.10514169375181734921e-12 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 31: {
      double t = 2*y100 - 63;
      return 0.30323425595617385705e0 + (0.95968346790597422934e-2 + (0.12931067776725883939e-3 + (-0.21938741702795543986e-5 + (-0.15202888584907373963e-6 + (0.61788350541116331411e-9 + (0.11957835742791248256e-9 + (-0.12598179834007710908e-11 - 0.75151817129574614194e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 32: {
      double t = 2*y100 - 65;
      return 0.32292521181517384379e0 + (0.10082957727001199408e-1 + (0.11257589426154962226e-3 + (-0.33670890319327881129e-5 + (-0.13910529040004008158e-6 + (0.19170714373047512945e-8 + (0.94840222377720494290e-10 + (-0.21650018351795353201e-11 - 0.37875211678024922689e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 33: {
      double t = 2*y100 - 67;
      return 0.34351233557911753862e0 + (0.10488575435572745309e-1 + (0.89209444197248726614e-4 + (-0.43893459576483345364e-5 + (-0.11488595830450424419e-6 + (0.28599494117122464806e-8 + (0.61537542799857777779e-10 - 0.24935749227658002212e-11 * t) * t) * t) * t) * t) * t) * t;
    }
    case 34: {
      double t = 2*y100 - 69;
      return 0.36480946642143669093e0 + (0.10789304203431861366e-1 + (0.60357993745283076834e-4 + (-0.51855862174130669389e-5 + (-0.83291664087289801313e-7 + (0.33898011178582671546e-8 + (0.27082948188277716482e-10 + (-0.23603379397408694974e-11 + 0.19328087692252869842e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 35: {
      double t = 2*y100 - 71;
      return 0.38658679935694939199e0 + (0.10966119158288804999e-1 + (0.27521612041849561426e-4 + (-0.57132774537670953638e-5 + (-0.48404772799207914899e-7 + (0.35268354132474570493e-8 + (-0.32383477652514618094e-11 + (-0.19334202915190442501e-11 + 0.32333189861286460270e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 36: {
      double t = 2*y100 - 73;
      return 0.40858275583808707870e0 + (0.11006378016848466550e-1 + (-0.76396376685213286033e-5 + (-0.59609835484245791439e-5 + (-0.13834610033859313213e-7 + (0.33406952974861448790e-8 + (-0.26474915974296612559e-10 + (-0.13750229270354351983e-11 + 0.36169366979417390637e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 37: {
      double t = 2*y100 - 75;
      return 0.43051714914006682977e0 + (0.10904106549500816155e-1 + (-0.43477527256787216909e-4 + (-0.59429739547798343948e-5 + (0.17639200194091885949e-7 + (0.29235991689639918688e-8 + (-0.41718791216277812879e-10 + (-0.81023337739508049606e-12 + 0.33618915934461994428e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 38: {
      double t = 2*y100 - 77;
      return 0.45210428135559607406e0 + (0.10659670756384400554e-1 + (-0.78488639913256978087e-4 + (-0.56919860886214735936e-5 + (0.44181850467477733407e-7 + (0.23694306174312688151e-8 + (-0.49492621596685443247e-10 + (-0.31827275712126287222e-12 + 0.27494438742721623654e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 39: {
      double t = 2*y100 - 79;
      return 0.47306491195005224077e0 + (0.10279006119745977570e-1 + (-0.11140268171830478306e-3 + (-0.52518035247451432069e-5 + (0.64846898158889479518e-7 + (0.17603624837787337662e-8 + (-0.51129481592926104316e-10 + (0.62674584974141049511e-13 + 0.20055478560829935356e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 40: {
      double t = 2*y100 - 81;
      return 0.49313638965719857647e0 + (0.97725799114772017662e-2 + (-0.14122854267291533334e-3 + (-0.46707252568834951907e-5 + (0.79421347979319449524e-7 + (0.11603027184324708643e-8 + (-0.48269605844397175946e-10 + (0.32477251431748571219e-12 + 0.12831052634143527985e-13 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 41: {
      double t = 2*y100 - 83;
      return 0.51208057433416004042e0 + (0.91542422354009224951e-2 + (-0.16726530230228647275e-3 + (-0.39964621752527649409e-5 + (0.88232252903213171454e-7 + (0.61343113364949928501e-9 + (-0.42516755603130443051e-10 + (0.47910437172240209262e-12 + 0.66784341874437478953e-14 * t) * t) * t) * t) * t) * t) * t) * t;
    }
    case 42: {
      double t = 2*y100 - 85;
      return 0.52968945458607484524e0 + (0.84400880445116786088e-2 + (-0.18908729783854258774e-3 + (-0.32725905467782951931e-5 + (0.91956190588652090659e-7 + (0.14593989152420122909e-9 + (-0.35239490687644444445e-10 + 0.54613829888448694898e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 43: {
      double t = 2*y100 - 87;
      return 0.54578857454330070965e0 + (0.76474155195880295311e-2 + (-0.20651230590808213884e-3 + (-0.25364339140543131706e-5 + (0.91455367999510681979e-7 + (-0.23061359005297528898e-9 + (-0.27512928625244444444e-10 + 0.54895806008493285579e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 44: {
      double t = 2*y100 - 89;
      return 0.56023851910298493910e0 + (0.67938321739997196804e-2 + (-0.21956066613331411760e-3 + (-0.18181127670443266395e-5 + (0.87650335075416845987e-7 + (-0.51548062050366615977e-9 + (-0.20068462174044444444e-10 + 0.50912654909758187264e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 45: {
      double t = 2*y100 - 91;
      return 0.57293478057455721150e0 + (0.58965321010394044087e-2 + (-0.22841145229276575597e-3 + (-0.11404605562013443659e-5 + (0.81430290992322326296e-7 + (-0.71512447242755357629e-9 + (-0.13372664928000000000e-10 + 0.44461498336689298148e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 46: {
      double t = 2*y100 - 93;
      return 0.58380635448407827360e0 + (0.49717469530842831182e-2 + (-0.23336001540009645365e-3 + (-0.51952064448608850822e-6 + (0.73596577815411080511e-7 + (-0.84020916763091566035e-9 + (-0.76700972702222222221e-11 + 0.36914462807972467044e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 47: {
      double t = 2*y100 - 95;
      return 0.59281340237769489597e0 + (0.40343592069379730568e-2 + (-0.23477963738658326185e-3 + (0.34615944987790224234e-7 + (0.64832803248395814574e-7 + (-0.90329163587627007971e-9 + (-0.30421940400000000000e-11 + 0.29237386653743536669e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 48: {
      double t = 2*y100 - 97;
      return 0.59994428743114271918e0 + (0.30976579788271744329e-2 + (-0.23308875765700082835e-3 + (0.51681681023846925160e-6 + (0.55694594264948268169e-7 + (-0.91719117313243464652e-9 + (0.53982743680000000000e-12 + 0.22050829296187771142e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 49: {
      double t = 2*y100 - 99;
      return 0.60521224471819875444e0 + (0.21732138012345456060e-2 + (-0.22872428969625997456e-3 + (0.92588959922653404233e-6 + (0.46612665806531930684e-7 + (-0.89393722514414153351e-9 + (0.31718550353777777778e-11 + 0.15705458816080549117e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 50: {
      double t = 2*y100 - 101;
      return 0.60865189969791123620e0 + (0.12708480848877451719e-2 + (-0.22212090111534847166e-3 + (0.12636236031532793467e-5 + (0.37904037100232937574e-7 + (-0.84417089968101223519e-9 + (0.49843180828444444445e-11 + 0.10355439441049048273e-12 * t) * t) * t) * t) * t) * t) * t;
    }
    case 51: {
      double t = 2*y100 - 103;
      return 0.61031580103499200191e0 + (0.39867436055861038223e-3 + (-0.21369573439579869291e-3 + (0.15339402129026183670e-5 + (0.29787479206646594442e-7 + (-0.77687792914228632974e-9 + (0.61192452741333333334e-11 + 0.60216691829459295780e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 52: {
      double t = 2*y100 - 105;
      return 0.61027109047879835868e0 + (-0.43680904508059878254e-3 + (-0.20383783788303894442e-3 + (0.17421743090883439959e-5 + (0.22400425572175715576e-7 + (-0.69934719320045128997e-9 + (0.67152759655111111110e-11 + 0.26419960042578359995e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 53: {
      double t = 2*y100 - 107;
      return 0.60859639489217430521e0 + (-0.12305921390962936873e-2 + (-0.19290150253894682629e-3 + (0.18944904654478310128e-5 + (0.15815530398618149110e-7 + (-0.61726850580964876070e-9 + 0.68987888999111111110e-11 * t) * t) * t) * t) * t) * t;
    }
    case 54: {
      double t = 2*y100 - 109;
      return 0.60537899426486075181e0 + (-0.19790062241395705751e-2 + (-0.18120271393047062253e-3 + (0.19974264162313241405e-5 + (0.10055795094298172492e-7 + (-0.53491997919318263593e-9 + (0.67794550295111111110e-11 - 0.17059208095741511603e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 55: {
      double t = 2*y100 - 111;
      return 0.60071229457904110537e0 + (-0.26795676776166354354e-2 + (-0.16901799553627508781e-3 + (0.20575498324332621581e-5 + (0.51077165074461745053e-8 + (-0.45536079828057221858e-9 + (0.64488005516444444445e-11 - 0.29311677573152766338e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 56: {
      double t = 2*y100 - 113;
      return 0.59469361520112714738e0 + (-0.33308208190600993470e-2 + (-0.15658501295912405679e-3 + (0.20812116912895417272e-5 + (0.93227468760614182021e-9 + (-0.38066673740116080415e-9 + (0.59806790359111111110e-11 - 0.36887077278950440597e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 57: {
      double t = 2*y100 - 115;
      return 0.58742228631775388268e0 + (-0.39321858196059227251e-2 + (-0.14410441141450122535e-3 + (0.20743790018404020716e-5 + (-0.25261903811221913762e-8 + (-0.31212416519526924318e-9 + (0.54328422462222222221e-11 - 0.40864152484979815972e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 58: {
      double t = 2*y100 - 117;
      return 0.57899804200033018447e0 + (-0.44838157005618913447e-2 + (-0.13174245966501437965e-3 + (0.20425306888294362674e-5 + (-0.53330296023875447782e-8 + (-0.25041289435539821014e-9 + (0.48490437205333333334e-11 - 0.42162206939169045177e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 59: {
      double t = 2*y100 - 119;
      return 0.56951968796931245974e0 + (-0.49864649488074868952e-2 + (-0.11963416583477567125e-3 + (0.19906021780991036425e-5 + (-0.75580140299436494248e-8 + (-0.19576060961919820491e-9 + (0.42613011928888888890e-11 - 0.41539443304115604377e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 60: {
      double t = 2*y100 - 121;
      return 0.55908401930063918964e0 + (-0.54413711036826877753e-2 + (-0.10788661102511914628e-3 + (0.19229663322982839331e-5 + (-0.92714731195118129616e-8 + (-0.14807038677197394186e-9 + (0.36920870298666666666e-11 - 0.39603726688419162617e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 61: {
      double t = 2*y100 - 123;
      return 0.54778496152925675315e0 + (-0.58501497933213396670e-2 + (-0.96582314317855227421e-4 + (0.18434405235069270228e-5 + (-0.10541580254317078711e-7 + (-0.10702303407788943498e-9 + (0.31563175582222222222e-11 - 0.36829748079110481422e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 62: {
      double t = 2*y100 - 125;
      return 0.53571290831682823999e0 + (-0.62147030670760791791e-2 + (-0.85782497917111760790e-4 + (0.17553116363443470478e-5 + (-0.11432547349815541084e-7 + (-0.72157091369041330520e-10 + (0.26630811607111111111e-11 - 0.33578660425893164084e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 63: {
      double t = 2*y100 - 127;
      return 0.52295422962048434978e0 + (-0.65371404367776320720e-2 + (-0.75530164941473343780e-4 + (0.16613725797181276790e-5 + (-0.12003521296598910761e-7 + (-0.42929753689181106171e-10 + (0.22170894940444444444e-11 - 0.30117697501065110505e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 64: {
      double t = 2*y100 - 129;
      return 0.50959092577577886140e0 + (-0.68197117603118591766e-2 + (-0.65852936198953623307e-4 + (0.15639654113906716939e-5 + (-0.12308007991056524902e-7 + (-0.18761997536910939570e-10 + (0.18198628922666666667e-11 - 0.26638355362285200932e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 65: {
      double t = 2*y100 - 131;
      return 0.49570040481823167970e0 + (-0.70647509397614398066e-2 + (-0.56765617728962588218e-4 + (0.14650274449141448497e-5 + (-0.12393681471984051132e-7 + (0.92904351801168955424e-12 + (0.14706755960177777778e-11 - 0.23272455351266325318e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 66: {
      double t = 2*y100 - 133;
      return 0.48135536250935238066e0 + (-0.72746293327402359783e-2 + (-0.48272489495730030780e-4 + (0.13661377309113939689e-5 + (-0.12302464447599382189e-7 + (0.16707760028737074907e-10 + (0.11672928324444444444e-11 - 0.20105801424709924499e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 67: {
      double t = 2*y100 - 135;
      return 0.46662374675511439448e0 + (-0.74517177649528487002e-2 + (-0.40369318744279128718e-4 + (0.12685621118898535407e-5 + (-0.12070791463315156250e-7 + (0.29105507892605823871e-10 + (0.90653314645333333334e-12 - 0.17189503312102982646e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 68: {
      double t = 2*y100 - 137;
      return 0.45156879030168268778e0 + (-0.75983560650033817497e-2 + (-0.33045110380705139759e-4 + (0.11732956732035040896e-5 + (-0.11729986947158201869e-7 + (0.38611905704166441308e-10 + (0.68468768305777777779e-12 - 0.14549134330396754575e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 69: {
      double t = 2*y100 - 139;
      return 0.43624909769330896904e0 + (-0.77168291040309554679e-2 + (-0.26283612321339907756e-4 + (0.10811018836893550820e-5 + (-0.11306707563739851552e-7 + (0.45670446788529607380e-10 + (0.49782492549333333334e-12 - 0.12191983967561779442e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 70: {
      double t = 2*y100 - 141;
      return 0.42071877443548481181e0 + (-0.78093484015052730097e-2 + (-0.20064596897224934705e-4 + (0.99254806680671890766e-6 + (-0.10823412088884741451e-7 + (0.50677203326904716247e-10 + (0.34200547594666666666e-12 - 0.10112698698356194618e-13 * t) * t) * t) * t) * t) * t) * t;
    }
    case 71: {
      double t = 2*y100 - 143;
      return 0.40502758809710844280e0 + (-0.78780384460872937555e-2 + (-0.14364940764532853112e-4 + (0.90803709228265217384e-6 + (-0.10298832847014466907e-7 + (0.53981671221969478551e-10 + (0.21342751381333333333e-12 - 0.82975901848387729274e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 72: {
      double t = 2*y100 - 145;
      return 0.38922115269731446690e0 + (-0.79249269708242064120e-2 + (-0.91595258799106970453e-5 + (0.82783535102217576495e-6 + (-0.97484311059617744437e-8 + (0.55889029041660225629e-10 + (0.10851981336888888889e-12 - 0.67278553237853459757e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 73: {
      double t = 2*y100 - 147;
      return 0.37334112915460307335e0 + (-0.79519385109223148791e-2 + (-0.44219833548840469752e-5 + (0.75209719038240314732e-6 + (-0.91848251458553190451e-8 + (0.56663266668051433844e-10 + (0.23995894257777777778e-13 - 0.53819475285389344313e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 74: {
      double t = 2*y100 - 149;
      return 0.35742543583374223085e0 + (-0.79608906571527956177e-2 + (-0.12530071050975781198e-6 + (0.68088605744900552505e-6 + (-0.86181844090844164075e-8 + (0.56530784203816176153e-10 + (-0.43120012248888888890e-13 - 0.42372603392496813810e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 75: {
      double t = 2*y100 - 151;
      return 0.34150846431979618536e0 + (-0.79534924968773806029e-2 + (0.37576885610891515813e-5 + (0.61419263633090524326e-6 + (-0.80565865409945960125e-8 + (0.55684175248749269411e-10 + (-0.95486860764444444445e-13 - 0.32712946432984510595e-14 * t) * t) * t) * t) * t) * t) * t;
    }
    case 76: {
      double t = 2*y100 - 153;
      return 0.32562129649136346824e0 + (-0.79313448067948884309e-2 + (0.72539159933545300034e-5 + (0.55195028297415503083e-6 + (-0.75063365335570475258e-8 + (0.54281686749699595941e-10 - 0.13545424295111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 77: {
      double t = 2*y100 - 155;
      return 0.30979191977078391864e0 + (-0.78959416264207333695e-2 + (0.10389774377677210794e-4 + (0.49404804463196316464e-6 + (-0.69722488229411164685e-8 + (0.52469254655951393842e-10 - 0.16507860650666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 78: {
      double t = 2*y100 - 157;
      return 0.29404543811214459904e0 + (-0.78486728990364155356e-2 + (0.13190885683106990459e-4 + (0.44034158861387909694e-6 + (-0.64578942561562616481e-8 + (0.50354306498006928984e-10 - 0.18614473550222222222e-12 * t) * t) * t) * t) * t) * t;
    }
    case 79: {
      double t = 2*y100 - 159;
      return 0.27840427686253660515e0 + (-0.77908279176252742013e-2 + (0.15681928798708548349e-4 + (0.39066226205099807573e-6 + (-0.59658144820660420814e-8 + (0.48030086420373141763e-10 - 0.20018995173333333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 80: {
      double t = 2*y100 - 161;
      return 0.26288838011163800908e0 + (-0.77235993576119469018e-2 + (0.17886516796198660969e-4 + (0.34482457073472497720e-6 + (-0.54977066551955420066e-8 + (0.45572749379147269213e-10 - 0.20852924954666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 81: {
      double t = 2*y100 - 163;
      return 0.24751539954181029717e0 + (-0.76480877165290370975e-2 + (0.19827114835033977049e-4 + (0.30263228619976332110e-6 + (-0.50545814570120129947e-8 + (0.43043879374212005966e-10 - 0.21228012028444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 82: {
      double t = 2*y100 - 165;
      return 0.23230087411688914593e0 + (-0.75653060136384041587e-2 + (0.21524991113020016415e-4 + (0.26388338542539382413e-6 + (-0.46368974069671446622e-8 + (0.40492715758206515307e-10 - 0.21238627815111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 83: {
      double t = 2*y100 - 167;
      return 0.21725840021297341931e0 + (-0.74761846305979730439e-2 + (0.23000194404129495243e-4 + (0.22837400135642906796e-6 + (-0.42446743058417541277e-8 + (0.37958104071765923728e-10 - 0.20963978568888888889e-12 * t) * t) * t) * t) * t) * t;
    }
    case 84: {
      double t = 2*y100 - 169;
      return 0.20239979200788191491e0 + (-0.73815761980493466516e-2 + (0.24271552727631854013e-4 + (0.19590154043390012843e-6 + (-0.38775884642456551753e-8 + (0.35470192372162901168e-10 - 0.20470131678222222222e-12 * t) * t) * t) * t) * t) * t;
    }
    case 85: {
      double t = 2*y100 - 171;
      return 0.18773523211558098962e0 + (-0.72822604530339834448e-2 + (0.25356688567841293697e-4 + (0.16626710297744290016e-6 + (-0.35350521468015310830e-8 + (0.33051896213898864306e-10 - 0.19811844544000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 86: {
      double t = 2*y100 - 173;
      return 0.17327341258479649442e0 + (-0.71789490089142761950e-2 + (0.26272046822383820476e-4 + (0.13927732375657362345e-6 + (-0.32162794266956859603e-8 + (0.30720156036105652035e-10 - 0.19034196304000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 87: {
      double t = 2*y100 - 175;
      return 0.15902166648328672043e0 + (-0.70722899934245504034e-2 + (0.27032932310132226025e-4 + (0.11474573347816568279e-6 + (-0.29203404091754665063e-8 + (0.28487010262547971859e-10 - 0.18174029063111111111e-12 * t) * t) * t) * t) * t) * t;
    }
    case 88: {
      double t = 2*y100 - 177;
      return 0.14498609036610283865e0 + (-0.69628725220045029273e-2 + (0.27653554229160596221e-4 + (0.92493727167393036470e-7 + (-0.26462055548683583849e-8 + (0.26360506250989943739e-10 - 0.17261211260444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 89: {
      double t = 2*y100 - 179;
      return 0.13117165798208050667e0 + (-0.68512309830281084723e-2 + (0.28147075431133863774e-4 + (0.72351212437979583441e-7 + (-0.23927816200314358570e-8 + (0.24345469651209833155e-10 - 0.16319736960000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 90: {
      double t = 2*y100 - 181;
      return 0.11758232561160626306e0 + (-0.67378491192463392927e-2 + (0.28525664781722907847e-4 + (0.54156999310046790024e-7 + (-0.21589405340123827823e-8 + (0.22444150951727334619e-10 - 0.15368675584000000000e-12 * t) * t) * t) * t) * t) * t;
    }
    case 91: {
      double t = 2*y100 - 183;
      return 0.10422112945361673560e0 + (-0.66231638959845581564e-2 + (0.28800551216363918088e-4 + (0.37758983397952149613e-7 + (-0.19435423557038933431e-8 + (0.20656766125421362458e-10 - 0.14422990012444444444e-12 * t) * t) * t) * t) * t) * t;
    }
    case 92: {
      double t = 2*y100 - 185;
      return 0.91090275493541084785e-1 + (-0.65075691516115160062e-2 + (0.28982078385527224867e-4 + (0.23014165807643012781e-7 + (-0.17454532910249875958e-8 + (0.18981946442680092373e-10 - 0.13494234691555555556e-12 * t) * t) * t) * t) * t) * t;
    }
    case 93: {
      double t = 2*y100 - 187;
      return 0.78191222288771379358e-1 + (-0.63914190297303976434e-2 + (0.29079759021299682675e-4 + (0.97885458059415717014e-8 + (-0.15635596116134296819e-8 + (0.17417110744051331974e-10 - 0.12591151763555555556e-12 * t) * t) * t) * t) * t) * t;
    }
    case 94: {
      double t = 2*y100 - 189;
      return 0.65524757106147402224e-1 + (-0.62750311956082444159e-2 + (0.29102328354323449795e-4 + (-0.20430838882727954582e-8 + (-0.13967781903855367270e-8 + (0.15958771833747057569e-10 - 0.11720175765333333333e-12 * t) * t) * t) * t) * t) * t;
    }
    case 95: {
      double t = 2*y100 - 191;
      return 0.53091065838453612773e-1 + (-0.61586898417077043662e-2 + (0.29057796072960100710e-4 + (-0.12597414620517987536e-7 + (-0.12440642607426861943e-8 + (0.14602787128447932137e-10 - 0.10885859114666666667e-12 * t) * t) * t) * t) * t) * t;
    }
    case 96: {
      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: 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.016931216931216931217))));
    }
  }
  /* Since 0 <= y100 < 101, this is only reached if x is NaN,
     in which case we should return NaN. */
  return NaN;
}

double FADDEEVA(w_im)(double x)
{
  if (x >= 0) {
    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;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
                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);
  }
  else { // = -FADDEEVA(w_im)(-x)
    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;
      /* 5-term expansion (rely on compiler for CSE), simplified from:
                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);
  }
}

/////////////////////////////////////////////////////////////////////////

// Compile with -DTEST_FADDEEVA to compile a little test program
#ifdef TEST_FADDEEVA

#ifdef __cplusplus
#  include <cstdio>
#else
#  include <stdio.h>
#endif

// compute relative error |b-a|/|a|, handling case of NaN and Inf,
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))
      return Inf; // "infinite" error
    return 0; // matching infinity/nan results counted as zero error
  }
  if (a == 0)
    return b == 0 ? 0 : Inf;
  else
    return fabs((b-a) / a);
}

int main(void) {
  double errmax_all = 0;
  {
    printf("############# w(z) tests #############\n");
#define NTST 57 // define instead of const for C compatibility
    cmplx z[NTST] = {
      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 */
      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) {
      cmplx fw = FADDEEVA(w)(z[i],0.);
      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);
      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("SUCCESS (max relative error = %g)\n", errmax);
    if (errmax > errmax_all) errmax_all = errmax;
  }
  {
#undef NTST
#define NTST 41 // define instead of const for C compatibility
    cmplx z[NTST] = {
      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
      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));                   \
      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)(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);
  }
  {
    // since erfi just calls through to erf, just one test should
    // 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] = { C(1.234,0.5678) };
    cmplx w[NTST] = { // erfi(z[i]), computed with Maple
      C(1.081032284405373149432716643834106923212,
        1.926775520840916645838949402886591180834)
    };
    TST(erfi, 0);
  }
  {
    // since erfcx just calls through to w, just one test should
    // 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] = { C(1.234,0.5678) };
    cmplx w[NTST] = { // erfcx(z[i]), computed with Maple
      C(0.3382187479799972294747793561190487832579,
        -0.1116077470811648467464927471872945833154)
    };
    TST(erfcx, 0);
  }
  {
#undef NTST
#define NTST 30 // define instead of const for C compatibility
    cmplx z[NTST] = {
      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
      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 48 // define instead of const for C compatibility
    cmplx z[NTST] = {
      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
      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);
  }
  printf("#####################################\n");
  printf("SUCCESS (max relative error = %g)\n", errmax_all);
}

#endif