GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Faddeeva.cc
Go to the documentation of this file.
1 // -*- mode:c++; tab-width:2; indent-tabs-mode:nil; -*-
2 
3 /* Copyright (c) 2012 Massachusetts Institute of Technology
4  *
5  * Permission is hereby granted, free of charge, to any person obtaining
6  * a copy of this software and associated documentation files (the
7  * "Software"), to deal in the Software without restriction, including
8  * without limitation the rights to use, copy, modify, merge, publish,
9  * distribute, sublicense, and/or sell copies of the Software, and to
10  * permit persons to whom the Software is furnished to do so, subject to
11  * the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be
14  * included in all copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20  * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21  * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22  * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23  */
24 
25 /* (Note that this file can be compiled with either C++, in which
26  case it uses C++ std::complex<double>, or C, in which case it
27  uses C99 double complex.) */
28 
29 /* Available at: http://ab-initio.mit.edu/Faddeeva
30 
31  Computes various error functions (erf, erfc, erfi, erfcx),
32  including the Dawson integral, in the complex plane, based
33  on algorithms for the computation of the Faddeeva function
34  w(z) = exp(-z^2) * erfc(-i*z).
35  Given w(z), the error functions are mostly straightforward
36  to compute, except for certain regions where we have to
37  switch to Taylor expansions to avoid cancellation errors
38  [e.g., near the origin for erf(z)].
39 
40  To compute the Faddeeva function, we use a combination of two
41  algorithms:
42 
43  For sufficiently large |z|, we use a continued-fraction expansion
44  for w(z) similar to those described in:
45 
46  Walter Gautschi, "Efficient computation of the complex error
47  function," SIAM J. Numer. Anal. 7(1), pp. 187-198 (1970)
48 
49  G. P. M. Poppe and C. M. J. Wijers, "More efficient computation
50  of the complex error function," ACM Trans. Math. Soft. 16(1),
51  pp. 38-46 (1990).
52 
53  Unlike those papers, however, we switch to a completely different
54  algorithm for smaller |z|:
55 
56  Mofreh R. Zaghloul and Ahmed N. Ali, "Algorithm 916: Computing the
57  Faddeyeva and Voigt Functions," ACM Trans. Math. Soft. 38(2), 15
58  (2011).
59 
60  (I initially used this algorithm for all z, but it turned out to be
61  significantly slower than the continued-fraction expansion for
62  larger |z|. On the other hand, it is competitive for smaller |z|,
63  and is significantly more accurate than the Poppe & Wijers code
64  in some regions, e.g., in the vicinity of z=1+1i.)
65 
66  Note that this is an INDEPENDENT RE-IMPLEMENTATION of these algorithms,
67  based on the description in the papers ONLY. In particular, I did
68  not refer to the authors' Fortran or Matlab implementations, respectively,
69  (which are under restrictive ACM copyright terms and therefore unusable
70  in free/open-source software).
71 
72  Steven G. Johnson, Massachusetts Institute of Technology
73  http://math.mit.edu/~stevenj
74  October 2012.
75 
76  -- Note that Algorithm 916 assumes that the erfc(x) function,
77  or rather the scaled function erfcx(x) = exp(x*x)*erfc(x),
78  is supplied for REAL arguments x. I originally used an
79  erfcx routine derived from DERFC in SLATEC, but I have
80  since replaced it with a much faster routine written by
81  me which uses a combination of continued-fraction expansions
82  and a lookup table of Chebyshev polynomials. For speed,
83  I implemented a similar algorithm for Im[w(x)] of real x,
84  since this comes up frequently in the other error functions.
85 
86  A small test program is included the end, which checks
87  the w(z) etc. results against several known values. To compile
88  the test function, compile with -DTEST_FADDEEVA (that is,
89  #define TEST_FADDEEVA).
90 
91  If HAVE_CONFIG_H is #defined (e.g., by compiling with -DHAVE_CONFIG_H),
92  then we #include "config.h", which is assumed to be a GNU autoconf-style
93  header defining HAVE_* macros to indicate the presence of features. In
94  particular, if HAVE_ISNAN and HAVE_ISINF are #defined, we use those
95  functions in math.h instead of defining our own, and if HAVE_ERF and/or
96  HAVE_ERFC are defined we use those functions from <cmath> for erf and
97  erfc of real arguments, respectively, instead of defining our own.
98 
99  REVISION HISTORY:
100  4 October 2012: Initial public release (SGJ)
101  5 October 2012: Revised (SGJ) to fix spelling error,
102  start summation for large x at round(x/a) (> 1)
103  rather than ceil(x/a) as in the original
104  paper, which should slightly improve performance
105  (and, apparently, slightly improves accuracy)
106  19 October 2012: Revised (SGJ) to fix bugs for large x, large -y,
107  and 15<x<26. Performance improvements. Prototype
108  now supplies default value for relerr.
109  24 October 2012: Switch to continued-fraction expansion for
110  sufficiently large z, for performance reasons.
111  Also, avoid spurious overflow for |z| > 1e154.
112  Set relerr argument to min(relerr,0.1).
113  27 October 2012: Enhance accuracy in Re[w(z)] taken by itself,
114  by switching to Alg. 916 in a region near
115  the real-z axis where continued fractions
116  have poor relative accuracy in Re[w(z)]. Thanks
117  to M. Zaghloul for the tip.
118  29 October 2012: Replace SLATEC-derived erfcx routine with
119  completely rewritten code by me, using a very
120  different algorithm which is much faster.
121  30 October 2012: Implemented special-case code for real z
122  (where real part is exp(-x^2) and imag part is
123  Dawson integral), using algorithm similar to erfx.
124  Export ImFaddeeva_w function to make Dawson's
125  integral directly accessible.
126  3 November 2012: Provide implementations of erf, erfc, erfcx,
127  and Dawson functions in Faddeeva:: namespace,
128  in addition to Faddeeva::w. Provide header
129  file Faddeeva.hh.
130  4 November 2012: Slightly faster erf for real arguments.
131  Updated MATLAB and Octave plugins.
132  27 November 2012: Support compilation with either C++ or
133  plain C (using C99 complex numbers).
134  For real x, use standard-library erf(x)
135  and erfc(x) if available (for C99 or C++11).
136  #include "config.h" if HAVE_CONFIG_H is #defined.
137  15 December 2012: Portability fixes (copysign, Inf/NaN creation),
138  use CMPLX/__builtin_complex if available in C,
139  slight accuracy improvements to erf and dawson
140  functions near the origin. Use gnulib functions
141  if GNULIB_NAMESPACE is defined.
142  18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson)
143 */
144 
145 /////////////////////////////////////////////////////////////////////////
146 /* If this file is compiled as a part of a larger project,
147  support using an autoconf-style config.h header file
148  (with various "HAVE_*" #defines to indicate features)
149  if HAVE_CONFIG_H is #defined (in GNU autotools style). */
150 
151 #if defined (HAVE_CONFIG_H)
152 # include "config.h"
153 #endif
154 
155 /////////////////////////////////////////////////////////////////////////
156 // macros to allow us to use either C++ or C (with C99 features)
157 
158 #if defined (__cplusplus)
159 
160 # include "lo-ieee.h"
161 
162 # include "Faddeeva.hh"
163 
164 # include <cfloat>
165 # include <cmath>
166 # include <limits>
167 
168 // use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS)
169 # define Inf octave::numeric_limits<double>::Inf ()
170 # define NaN octave::numeric_limits<double>::NaN ()
171 
172 typedef std::complex<double> cmplx;
173 
174 // Use C-like complex syntax, since the C syntax is more restrictive
175 # define cexp(z) exp(z)
176 # define creal(z) real(z)
177 # define cimag(z) imag(z)
178 # define cpolar(r,t) polar(r,t)
179 
180 # define C(a,b) cmplx(a,b)
181 
182 # define FADDEEVA(name) Faddeeva::name
183 # define FADDEEVA_RE(name) Faddeeva::name
184 
185 // isnan/isinf were introduced in C++11
186 # if defined (lo_ieee_isnan) && defined (lo_ieee_isinf)
187 # define isnan lo_ieee_isnan
188 # define isinf lo_ieee_isinf
189 # elif (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF))
190 static inline bool my_isnan(double x) { return x != x; }
191 # define isnan my_isnan
192 static inline bool my_isinf(double x) { return 1/x == 0.; }
193 # define isinf my_isinf
194 # elif (__cplusplus >= 201103L)
195 // g++ gets confused between the C and C++ isnan/isinf functions
196 # define isnan std::isnan
197 # define isinf std::isinf
198 # endif
199 
200 // copysign was introduced in C++11 (and is also in POSIX and C99)
201 # if defined(_WIN32) || defined(__WIN32__)
202 # define copysign _copysign // of course MS had to be different
203 # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX)
204 static inline double my_copysign(double x, double y) { return y<0 ? -x : x; }
205 # define copysign my_copysign
206 # endif
207 
208 #else // !__cplusplus, i.e., pure C (requires C99 features)
209 
210 # include "Faddeeva.h"
211 
212 # define _GNU_SOURCE // enable GNU libc NAN extension if possible
213 
214 # include <float.h>
215 # include <math.h>
216 
217 typedef double complex cmplx;
218 
219 # define FADDEEVA(name) Faddeeva_ ## name
220 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
221 
222 /* Constructing complex numbers like 0+i*NaN is problematic in C99
223  without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
224  I is a complex (rather than imaginary) constant. For some reason,
225  however, it works fine in (pre-4.7) gcc if I define Inf and NaN as
226  1/0 and 0/0 (and only if I compile with optimization -O1 or more),
227  but not if I use the INFINITY or NAN macros. */
228 
229 /* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro
230  may not be defined unless we are using a recent (2012) version of
231  glibc and compile with -std=c11... note that icc lies about being
232  gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
233 # if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER))
234 # define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
235 # endif
236 
237 # if defined (CMPLX) // C11
238 # define C(a,b) CMPLX(a,b)
239 # define Inf INFINITY // C99 infinity
240 # if defined (NAN) // GNU libc extension
241 # define NaN NAN
242 # else
243 # define NaN (0./0.) // NaN
244 # endif
245 # else
246 # define C(a,b) ((a) + I*(b))
247 # define Inf (1./0.)
248 # define NaN (0./0.)
249 # endif
250 
251 static inline cmplx cpolar(double r, double t)
252 {
253  if (r == 0.0 && !isnan(t))
254  return 0.0;
255  else
256  return C(r * cos(t), r * sin(t));
257 }
258 
259 #endif // !__cplusplus, i.e., pure C (requires C99 features)
260 
261 /////////////////////////////////////////////////////////////////////////
262 // Auxiliary routines to compute other special functions based on w(z)
263 
264 // compute erfcx(z) = exp(z^2) erfz(z)
265 cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
266 {
267  return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
268 }
269 
270 // compute the error function erf(x)
271 double FADDEEVA_RE(erf)(double x)
272 {
273 #if !defined(__cplusplus)
274  return erf(x); // C99 supplies erf in math.h
275 #elif (__cplusplus >= 201103L) || defined(HAVE_ERF)
276  return ::erf(x); // C++11 supplies std::erf in cmath
277 #else
278  double mx2 = -x*x;
279  if (mx2 < -750) // underflow
280  return (x >= 0 ? 1.0 : -1.0);
281 
282  if (x >= 0) {
283  if (x < 8e-2) goto taylor;
284  return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
285  }
286  else { // x < 0
287  if (x > -8e-2) goto taylor;
288  return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0;
289  }
290 
291  // Use Taylor series for small |x|, to avoid cancellation inaccuracy
292  // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...)
293  taylor:
294  return x * (1.1283791670955125739
295  + mx2 * (0.37612638903183752464
296  + mx2 * (0.11283791670955125739
297  + mx2 * (0.026866170645131251760
298  + mx2 * 0.0052239776254421878422))));
299 #endif
300 }
301 
302 // compute the error function erf(z)
303 cmplx FADDEEVA(erf)(cmplx z, double relerr)
304 {
305  double x = creal(z), y = cimag(z);
306 
307  if (y == 0)
308  return C(FADDEEVA_RE(erf)(x),
309  y); // preserve sign of 0
310  if (x == 0) // handle separately for speed & handling of y = Inf or NaN
311  return C(x, // preserve sign of 0
312  /* handle y -> Inf limit manually, since
313  exp(y^2) -> Inf but Im[w(y)] -> 0, so
314  IEEE will give us a NaN when it should be Inf */
315  y*y > 720 ? (y > 0 ? Inf : -Inf)
316  : exp(y*y) * FADDEEVA(w_im)(y));
317 
318  double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
319  double mIm_z2 = -2*x*y; // Im(-z^2)
320  if (mRe_z2 < -750) // underflow
321  return (x >= 0 ? 1.0 : -1.0);
322 
323  /* Handle positive and negative x via different formulas,
324  using the mirror symmetries of w, to avoid overflow/underflow
325  problems from multiplying exponentially large and small quantities. */
326  if (x >= 0) {
327  if (x < 8e-2) {
328  if (fabs(y) < 1e-2)
329  goto taylor;
330  else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
331  goto taylor_erfi;
332  }
333  /* don't use complex exp function, since that will produce spurious NaN
334  values when multiplying w in an overflow situation. */
335  return 1.0 - exp(mRe_z2) *
336  (C(cos(mIm_z2), sin(mIm_z2))
337  * FADDEEVA(w)(C(-y,x), relerr));
338  }
339  else { // x < 0
340  if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
341  if (fabs(y) < 1e-2)
342  goto taylor;
343  else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
344  goto taylor_erfi;
345  }
346  else if (isnan(x))
347  return C(NaN, y == 0 ? 0 : NaN);
348  /* don't use complex exp function, since that will produce spurious NaN
349  values when multiplying w in an overflow situation. */
350  return exp(mRe_z2) *
351  (C(cos(mIm_z2), sin(mIm_z2))
352  * FADDEEVA(w)(C(y,-x), relerr)) - 1.0;
353  }
354 
355  // Use Taylor series for small |z|, to avoid cancellation inaccuracy
356  // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
357  taylor:
358  {
359  cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
360  return z * (1.1283791670955125739
361  + mz2 * (0.37612638903183752464
362  + mz2 * (0.11283791670955125739
363  + mz2 * (0.026866170645131251760
364  + mz2 * 0.0052239776254421878422))));
365  }
366 
367  /* for small |x| and small |xy|,
368  use Taylor series to avoid cancellation inaccuracy:
369  erf(x+iy) = erf(iy)
370  + 2*exp(y^2)/sqrt(pi) *
371  [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
372  - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
373  where:
374  erf(iy) = exp(y^2) * Im[w(y)]
375  */
376  taylor_erfi:
377  {
378  double x2 = x*x, y2 = y*y;
379  double expy2 = exp(y2);
380  return C
381  (expy2 * x * (1.1283791670955125739
382  - x2 * (0.37612638903183752464
383  + 0.75225277806367504925*y2)
384  + x2*x2 * (0.11283791670955125739
385  + y2 * (0.45135166683820502956
386  + 0.15045055561273500986*y2))),
387  expy2 * (FADDEEVA(w_im)(y)
388  - x2*y * (1.1283791670955125739
389  - x2 * (0.56418958354775628695
390  + 0.37612638903183752464*y2))));
391  }
392 }
393 
394 // erfi(z) = -i erf(iz)
395 cmplx FADDEEVA(erfi)(cmplx z, double relerr)
396 {
397  cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr);
398  return C(cimag(e), -creal(e));
399 }
400 
401 // erfi(x) = -i erf(ix)
402 double FADDEEVA_RE(erfi)(double x)
403 {
404  return x*x > 720 ? (x > 0 ? Inf : -Inf)
405  : exp(x*x) * FADDEEVA(w_im)(x);
406 }
407 
408 // erfc(x) = 1 - erf(x)
409 double FADDEEVA_RE(erfc)(double x)
410 {
411 #if !defined(__cplusplus)
412  return erfc(x); // C99 supplies erfc in math.h
413 #elif (__cplusplus >= 201103L) || defined(HAVE_ERFC)
414  return ::erfc(x); // C++11 supplies std::erfc in cmath
415 #else
416  if (x*x > 750) // underflow
417  return (x >= 0 ? 0.0 : 2.0);
418  return x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x)
419  : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x);
420 #endif
421 }
422 
423 // erfc(z) = 1 - erf(z)
424 cmplx FADDEEVA(erfc)(cmplx z, double relerr)
425 {
426  double x = creal(z), y = cimag(z);
427 
428  if (x == 0.)
429  return C(1,
430  /* handle y -> Inf limit manually, since
431  exp(y^2) -> Inf but Im[w(y)] -> 0, so
432  IEEE will give us a NaN when it should be Inf */
433  y*y > 720 ? (y > 0 ? -Inf : Inf)
434  : -exp(y*y) * FADDEEVA(w_im)(y));
435  if (y == 0.) {
436  if (x*x > 750) // underflow
437  return C(x >= 0 ? 0.0 : 2.0,
438  -y); // preserve sign of 0
439  return C(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x)
440  : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x),
441  -y); // preserve sign of zero
442  }
443 
444  double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
445  double mIm_z2 = -2*x*y; // Im(-z^2)
446  if (mRe_z2 < -750) // underflow
447  return (x >= 0 ? 0.0 : 2.0);
448 
449  if (x >= 0)
450  return cexp(C(mRe_z2, mIm_z2))
451  * FADDEEVA(w)(C(-y,x), relerr);
452  else
453  return 2.0 - cexp(C(mRe_z2, mIm_z2))
454  * FADDEEVA(w)(C(y,-x), relerr);
455 }
456 
457 // compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x)
458 double FADDEEVA_RE(Dawson)(double x)
459 {
460  const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
461  return spi2 * FADDEEVA(w_im)(x);
462 }
463 
464 // compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)
465 cmplx FADDEEVA(Dawson)(cmplx z, double relerr)
466 {
467  const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
468  double x = creal(z), y = cimag(z);
469 
470  // handle axes separately for speed & proper handling of x or y = Inf or NaN
471  if (y == 0)
472  return C(spi2 * FADDEEVA(w_im)(x),
473  -y); // preserve sign of 0
474  if (x == 0) {
475  double y2 = y*y;
476  if (y2 < 2.5e-5) { // Taylor expansion
477  return C(x, // preserve sign of 0
478  y * (1.
479  + y2 * (0.6666666666666666666666666666666666666667
480  + y2 * 0.26666666666666666666666666666666666667)));
481  }
482  return C(x, // preserve sign of 0
483  spi2 * (y >= 0
484  ? exp(y2) - FADDEEVA_RE(erfcx)(y)
485  : FADDEEVA_RE(erfcx)(-y) - exp(y2)));
486  }
487 
488  double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
489  double mIm_z2 = -2*x*y; // Im(-z^2)
490  cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
491 
492  /* Handle positive and negative x via different formulas,
493  using the mirror symmetries of w, to avoid overflow/underflow
494  problems from multiplying exponentially large and small quantities. */
495  if (y >= 0) {
496  if (y < 5e-3) {
497  if (fabs(x) < 5e-3)
498  goto taylor;
499  else if (fabs(mIm_z2) < 5e-3)
500  goto taylor_realaxis;
501  }
502  cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr);
503  return spi2 * C(-cimag(res), creal(res));
504  }
505  else { // y < 0
506  if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
507  if (fabs(x) < 5e-3)
508  goto taylor;
509  else if (fabs(mIm_z2) < 5e-3)
510  goto taylor_realaxis;
511  }
512  else if (isnan(y))
513  return C(x == 0 ? 0 : NaN, NaN);
514  cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2);
515  return spi2 * C(-cimag(res), creal(res));
516  }
517 
518  // Use Taylor series for small |z|, to avoid cancellation inaccuracy
519  // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
520  taylor:
521  return z * (1.
522  + mz2 * (0.6666666666666666666666666666666666666667
523  + mz2 * 0.2666666666666666666666666666666666666667));
524 
525  /* for small |y| and small |xy|,
526  use Taylor series to avoid cancellation inaccuracy:
527  dawson(x + iy)
528  = D + y^2 (D + x - 2Dx^2)
529  + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
530  + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
531  + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
532  where D = dawson(x)
533 
534  However, for large |x|, 2Dx -> 1 which gives cancellation problems in
535  this series (many of the leading terms cancel). So, for large |x|,
536  we need to substitute a continued-fraction expansion for D.
537 
538  dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))
539 
540  The 6 terms shown here seems to be the minimum needed to be
541  accurate as soon as the simpler Taylor expansion above starts
542  breaking down. Using this 6-term expansion, factoring out the
543  denominator, and simplifying with Maple, we obtain:
544 
545  Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x
546  = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
547  Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y
548  = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4
549 
550  Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction
551  expansion for the real part, and a 2-term expansion for the imaginary
552  part. (This avoids overflow problems for huge |x|.) This yields:
553 
554  Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
555  Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)
556 
557  */
558  taylor_realaxis:
559  {
560  double x2 = x*x;
561  if (x2 > 1600) { // |x| > 40
562  double y2 = y*y;
563  if (x2 > 25e14) {// |x| > 5e7
564  double xy2 = (x*y)*(x*y);
565  return C((0.5 + y2 * (0.5 + 0.25*y2
566  - 0.16666666666666666667*xy2)) / x,
567  y * (-1 + y2 * (-0.66666666666666666667
568  + 0.13333333333333333333*xy2
569  - 0.26666666666666666667*y2))
570  / (2*x2 - 1));
571  }
572  return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
573  C(x * (33 + x2 * (-28 + 4*x2)
574  + y2 * (18 - 4*x2 + 4*y2)),
575  y * (-15 + x2 * (24 - 4*x2)
576  + y2 * (4*x2 - 10 - 4*y2)));
577  }
578  else {
579  double D = spi2 * FADDEEVA(w_im)(x);
580  double y2 = y*y;
581  return C
582  (D + y2 * (D + x - 2*D*x2)
583  + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
584  + x * (0.83333333333333333333
585  - 0.33333333333333333333 * x2)),
586  y * (1 - 2*D*x
587  + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
588  + y2*y2 * (0.26666666666666666667 -
589  x2 * (0.6 - 0.13333333333333333333 * x2)
590  - D*x * (1 - x2 * (1.3333333333333333333
591  - 0.26666666666666666667 * x2)))));
592  }
593  }
594 }
595 
596 /////////////////////////////////////////////////////////////////////////
597 
598 // return sinc(x) = sin(x)/x, given both x and sin(x)
599 // [since we only use this in cases where sin(x) has already been computed]
600 static inline double sinc(double x, double sinx) {
601  return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x;
602 }
603 
604 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2
605 static inline double sinh_taylor(double x) {
606  return x * (1 + (x*x) * (0.1666666666666666666667
607  + 0.00833333333333333333333 * (x*x)));
608 }
609 
610 static inline double sqr(double x) { return x*x; }
611 
612 // precomputed table of expa2n2[n-1] = exp(-a2*n*n)
613 // for double-precision a2 = 0.26865... in FADDEEVA(w), below.
614 static const double expa2n2[] = {
615  7.64405281671221563e-01,
616  3.41424527166548425e-01,
617  8.91072646929412548e-02,
618  1.35887299055460086e-02,
619  1.21085455253437481e-03,
620  6.30452613933449404e-05,
621  1.91805156577114683e-06,
622  3.40969447714832381e-08,
623  3.54175089099469393e-10,
624  2.14965079583260682e-12,
625  7.62368911833724354e-15,
626  1.57982797110681093e-17,
627  1.91294189103582677e-20,
628  1.35344656764205340e-23,
629  5.59535712428588720e-27,
630  1.35164257972401769e-30,
631  1.90784582843501167e-34,
632  1.57351920291442930e-38,
633  7.58312432328032845e-43,
634  2.13536275438697082e-47,
635  3.51352063787195769e-52,
636  3.37800830266396920e-57,
637  1.89769439468301000e-62,
638  6.22929926072668851e-68,
639  1.19481172006938722e-73,
640  1.33908181133005953e-79,
641  8.76924303483223939e-86,
642  3.35555576166254986e-92,
643  7.50264110688173024e-99,
644  9.80192200745410268e-106,
645  7.48265412822268959e-113,
646  3.33770122566809425e-120,
647  8.69934598159861140e-128,
648  1.32486951484088852e-135,
649  1.17898144201315253e-143,
650  6.13039120236180012e-152,
651  1.86258785950822098e-160,
652  3.30668408201432783e-169,
653  3.43017280887946235e-178,
654  2.07915397775808219e-187,
655  7.36384545323984966e-197,
656  1.52394760394085741e-206,
657  1.84281935046532100e-216,
658  1.30209553802992923e-226,
659  5.37588903521080531e-237,
660  1.29689584599763145e-247,
661  1.82813078022866562e-258,
662  1.50576355348684241e-269,
663  7.24692320799294194e-281,
664  2.03797051314726829e-292,
665  3.34880215927873807e-304,
666  0.0 // underflow (also prevents reads past array end, below)
667 };
668 
669 /////////////////////////////////////////////////////////////////////////
670 
671 cmplx FADDEEVA(w)(cmplx z, double relerr)
672 {
673  if (creal(z) == 0.0)
674  return C(FADDEEVA_RE(erfcx)(cimag(z)),
675  creal(z)); // give correct sign of 0 in cimag(w)
676  else if (cimag(z) == 0)
677  return C(exp(-sqr(creal(z))),
678  FADDEEVA(w_im)(creal(z)));
679 
680  double a, a2, c;
681  if (relerr <= DBL_EPSILON) {
682  relerr = DBL_EPSILON;
683  a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))
684  c = 0.329973702884629072537; // (2/pi) * a;
685  a2 = 0.268657157075235951582; // a^2
686  }
687  else {
688  const double pi = 3.14159265358979323846264338327950288419716939937510582;
689  if (relerr > 0.1) relerr = 0.1; // not sensible to compute < 1 digit
690  a = pi / sqrt(-log(relerr*0.5));
691  c = (2/pi)*a;
692  a2 = a*a;
693  }
694  const double x = fabs(creal(z));
695  const double y = cimag(z), ya = fabs(y);
696 
697  cmplx ret = 0.; // return value
698 
699  double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;
700 
701 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z|
702 
703 #if USE_CONTINUED_FRACTION
704  if (ya > 7 || (x > 6 // continued fraction is faster
705  /* As pointed out by M. Zaghloul, the continued
706  fraction seems to give a large relative error in
707  Re w(z) for |x| ~ 6 and small |y|, so use
708  algorithm 816 in this region: */
709  && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
710 
711  /* Poppe & Wijers suggest using a number of terms
712  nu = 3 + 1442 / (26*rho + 77)
713  where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
714  (They only use this expansion for rho >= 1, but rho a little less
715  than 1 seems okay too.)
716  Instead, I did my own fit to a slightly different function
717  that avoids the hypotenuse calculation, using NLopt to minimize
718  the sum of the squares of the errors in nu with the constraint
719  that the estimated nu be >= minimum nu to attain machine precision.
720  I also separate the regions where nu == 2 and nu == 1. */
721  const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
722  double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
723  if (x + ya > 4000) { // nu <= 2
724  if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z
725  // scale to avoid overflow
726  if (x > ya) {
727  double yax = ya / xs;
728  double denom = ispi / (xs + yax*ya);
729  ret = C(denom*yax, denom);
730  }
731  else if (isinf(ya))
732  return ((isnan(x) || y < 0)
733  ? C(NaN,NaN) : C(0,0));
734  else {
735  double xya = xs / ya;
736  double denom = ispi / (xya*xs + ya);
737  ret = C(denom, denom*xya);
738  }
739  }
740  else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
741  double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
742  double denom = ispi / (dr*dr + di*di);
743  ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
744  }
745  }
746  else { // compute nu(z) estimate and do general continued fraction
747  const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit
748  double nu = floor(c0 + c1 / (c2*x + c3*ya + c4));
749  double wr = xs, wi = ya;
750  for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) {
751  // w <- z - nu/w:
752  double denom = nu / (wr*wr + wi*wi);
753  wr = xs - wr * denom;
754  wi = ya + wi * denom;
755  }
756  { // w(z) = i/sqrt(pi) / w:
757  double denom = ispi / (wr*wr + wi*wi);
758  ret = C(denom*wi, denom*wr);
759  }
760  }
761  if (y < 0) {
762  // use w(z) = 2.0*exp(-z*z) - w(-z),
763  // but be careful of overflow in exp(-z*z)
764  // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
765  return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
766  }
767  else
768  return ret;
769  }
770 #else // !USE_CONTINUED_FRACTION
771  if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision
772  const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
773  double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
774  // scale to avoid overflow
775  if (x > ya) {
776  double yax = ya / xs;
777  double denom = ispi / (xs + yax*ya);
778  ret = C(denom*yax, denom);
779  }
780  else {
781  double xya = xs / ya;
782  double denom = ispi / (xya*xs + ya);
783  ret = C(denom, denom*xya);
784  }
785  if (y < 0) {
786  // use w(z) = 2.0*exp(-z*z) - w(-z),
787  // but be careful of overflow in exp(-z*z)
788  // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
789  return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
790  }
791  else
792  return ret;
793  }
794 #endif // !USE_CONTINUED_FRACTION
795 
796  /* Note: The test that seems to be suggested in the paper is x <
797  sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2)
798  underflows to zero and sum1,sum2,sum4 are zero. However, long
799  before this occurs, the sum1,sum2,sum4 contributions are
800  negligible in double precision; I find that this happens for x >
801  about 6, for all y. On the other hand, I find that the case
802  where we compute all of the sums is faster (at least with the
803  precomputed expa2n2 table) until about x=10. Furthermore, if we
804  try to compute all of the sums for x > 20, I find that we
805  sometimes run into numerical problems because underflow/overflow
806  problems start to appear in the various coefficients of the sums,
807  below. Therefore, we use x < 10 here. */
808  else if (x < 10) {
809  double prod2ax = 1, prodm2ax = 1;
810  double expx2;
811 
812  if (isnan(y))
813  return C(y,y);
814 
815  /* Somewhat ugly copy-and-paste duplication here, but I see significant
816  speedups from using the special-case code with the precomputed
817  exponential, and the x < 5e-4 special case is needed for accuracy. */
818 
819  if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table
820  if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
821  const double x2 = x*x;
822  expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
823  // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
824  const double ax2 = 1.036642960860171859744*x; // 2*a*x
825  const double exp2ax =
826  1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
827  const double expm2ax =
828  1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
829  for (int n = 1; 1; ++n) {
830  const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
831  prod2ax *= exp2ax;
832  prodm2ax *= expm2ax;
833  sum1 += coef;
834  sum2 += coef * prodm2ax;
835  sum3 += coef * prod2ax;
836 
837  // really = sum5 - sum4
838  sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
839 
840  // test convergence via sum3
841  if (coef * prod2ax < relerr * sum3) break;
842  }
843  }
844  else { // x > 5e-4, compute sum4 and sum5 separately
845  expx2 = exp(-x*x);
846  const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
847  for (int n = 1; 1; ++n) {
848  const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
849  prod2ax *= exp2ax;
850  prodm2ax *= expm2ax;
851  sum1 += coef;
852  sum2 += coef * prodm2ax;
853  sum4 += (coef * prodm2ax) * (a*n);
854  sum3 += coef * prod2ax;
855  sum5 += (coef * prod2ax) * (a*n);
856  // test convergence via sum5, since this sum has the slowest decay
857  if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
858  }
859  }
860  }
861  else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly
862  const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
863  if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
864  const double x2 = x*x;
865  expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
866  for (int n = 1; 1; ++n) {
867  const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
868  prod2ax *= exp2ax;
869  prodm2ax *= expm2ax;
870  sum1 += coef;
871  sum2 += coef * prodm2ax;
872  sum3 += coef * prod2ax;
873 
874  // really = sum5 - sum4
875  sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
876 
877  // test convergence via sum3
878  if (coef * prod2ax < relerr * sum3) break;
879  }
880  }
881  else { // x > 5e-4, compute sum4 and sum5 separately
882  expx2 = exp(-x*x);
883  for (int n = 1; 1; ++n) {
884  const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
885  prod2ax *= exp2ax;
886  prodm2ax *= expm2ax;
887  sum1 += coef;
888  sum2 += coef * prodm2ax;
889  sum4 += (coef * prodm2ax) * (a*n);
890  sum3 += coef * prod2ax;
891  sum5 += (coef * prod2ax) * (a*n);
892  // test convergence via sum5, since this sum has the slowest decay
893  if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
894  }
895  }
896  }
897  const double expx2erfcxy = // avoid spurious overflow for large negative y
898  y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision
899  ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x);
900  if (y > 5) { // imaginary terms cancel
901  const double sinxy = sin(x*y);
902  ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
903  + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
904  }
905  else {
906  double xs = creal(z);
907  const double sinxy = sin(xs*y);
908  const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y);
909  const double coef1 = expx2erfcxy - c*y*sum1;
910  const double coef2 = c*xs*expx2;
911  ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
912  coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
913  }
914  }
915  else { // x large: only sum3 & sum5 contribute (see above note)
916  if (isnan(x))
917  return C(x,x);
918  if (isnan(y))
919  return C(y,y);
920 
921 #if USE_CONTINUED_FRACTION
922  ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
923 #else
924  if (y < 0) {
925  /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so
926  erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
927  if y*y - x*x > -36 or so. So, compute this term just in case.
928  We also need the -exp(-x*x) term to compute Re[w] accurately
929  in the case where y is very small. */
930  ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
931  }
932  else
933  ret = exp(-x*x); // not negligible in real part if y very small
934 #endif
935  // (round instead of ceil as in original paper; note that x/a > 1 here)
936  double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0
937  double dx = a*n0 - x;
938  sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y);
939  sum5 = a*n0 * sum3;
940  double exp1 = exp(4*a*dx), exp1dn = 1;
941  int dn;
942  for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms
943  double np = n0 + dn, nm = n0 - dn;
944  double tp = exp(-sqr(a*dn+dx));
945  double tm = tp * (exp1dn *= exp1); // trick to get tm from tp
946  tp /= (a2*(np*np) + y*y);
947  tm /= (a2*(nm*nm) + y*y);
948  sum3 += tp + tm;
949  sum5 += a * (np * tp + nm * tm);
950  if (a * (np * tp + nm * tm) < relerr * sum5) goto finish;
951  }
952  while (1) { // loop over n0+dn terms only (since n0-dn <= 0)
953  double np = n0 + dn++;
954  double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y);
955  sum3 += tp;
956  sum5 += a * np * tp;
957  if (a * np * tp < relerr * sum5) goto finish;
958  }
959  }
960  finish:
961  return ret + C((0.5*c)*y*(sum2+sum3),
962  (0.5*c)*copysign(sum5-sum4, creal(z)));
963 }
964 
965 /////////////////////////////////////////////////////////////////////////
966 
967 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by
968  Steven G. Johnson, October 2012.
969 
970  This function combines a few different ideas.
971 
972  First, for x > 50, it uses a continued-fraction expansion (same as
973  for the Faddeeva function, but with algebraic simplifications for z=i*x).
974 
975  Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations,
976  but with two twists:
977 
978  a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation,
979  inspired by a similar transformation in the octave-forge/specfun
980  erfcx by Soren Hauberg, results in much faster Chebyshev convergence
981  than other simple transformations I have examined.
982 
983  b) Instead of using a single Chebyshev polynomial for the entire
984  [0,1] y interval, we break the interval up into 100 equal
985  subintervals, with a switch/lookup table, and use much lower
986  degree Chebyshev polynomials in each subinterval. This greatly
987  improves performance in my tests.
988 
989  For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
990  with the usual checks for overflow etcetera.
991 
992  Performance-wise, it seems to be substantially faster than either
993  the SLATEC DERFC function [or an erfcx function derived therefrom]
994  or Cody's CALERF function (from netlib.org/specfun), while
995  retaining near machine precision in accuracy. */
996 
997 /* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x).
998 
999  Uses a look-up table of 100 different Chebyshev polynomials
1000  for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
1001  with the help of Maple and a little shell script. This allows
1002  the Chebyshev polynomials to be of significantly lower degree (about 1/4)
1003  compared to fitting the whole [0,1] interval with a single polynomial. */
1004 static double erfcx_y100(double y100)
1005 {
1006  switch (static_cast<int> (y100)) {
1007 case 0: {
1008 double t = 2*y100 - 1;
1009 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;
1010 }
1011 case 1: {
1012 double t = 2*y100 - 3;
1013 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;
1014 }
1015 case 2: {
1016 double t = 2*y100 - 5;
1017 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;
1018 }
1019 case 3: {
1020 double t = 2*y100 - 7;
1021 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;
1022 }
1023 case 4: {
1024 double t = 2*y100 - 9;
1025 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;
1026 }
1027 case 5: {
1028 double t = 2*y100 - 11;
1029 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;
1030 }
1031 case 6: {
1032 double t = 2*y100 - 13;
1033 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;
1034 }
1035 case 7: {
1036 double t = 2*y100 - 15;
1037 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;
1038 }
1039 case 8: {
1040 double t = 2*y100 - 17;
1041 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;
1042 }
1043 case 9: {
1044 double t = 2*y100 - 19;
1045 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;
1046 }
1047 case 10: {
1048 double t = 2*y100 - 21;
1049 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;
1050 }
1051 case 11: {
1052 double t = 2*y100 - 23;
1053 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;
1054 }
1055 case 12: {
1056 double t = 2*y100 - 25;
1057 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;
1058 }
1059 case 13: {
1060 double t = 2*y100 - 27;
1061 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;
1062 }
1063 case 14: {
1064 double t = 2*y100 - 29;
1065 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;
1066 }
1067 case 15: {
1068 double t = 2*y100 - 31;
1069 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;
1070 }
1071 case 16: {
1072 double t = 2*y100 - 33;
1073 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;
1074 }
1075 case 17: {
1076 double t = 2*y100 - 35;
1077 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;
1078 }
1079 case 18: {
1080 double t = 2*y100 - 37;
1081 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;
1082 }
1083 case 19: {
1084 double t = 2*y100 - 39;
1085 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;
1086 }
1087 case 20: {
1088 double t = 2*y100 - 41;
1089 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;
1090 }
1091 case 21: {
1092 double t = 2*y100 - 43;
1093 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;
1094 }
1095 case 22: {
1096 double t = 2*y100 - 45;
1097 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;
1098 }
1099 case 23: {
1100 double t = 2*y100 - 47;
1101 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;
1102 }
1103 case 24: {
1104 double t = 2*y100 - 49;
1105 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;
1106 }
1107 case 25: {
1108 double t = 2*y100 - 51;
1109 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;
1110 }
1111 case 26: {
1112 double t = 2*y100 - 53;
1113 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;
1114 }
1115 case 27: {
1116 double t = 2*y100 - 55;
1117 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;
1118 }
1119 case 28: {
1120 double t = 2*y100 - 57;
1121 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;
1122 }
1123 case 29: {
1124 double t = 2*y100 - 59;
1125 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;
1126 }
1127 case 30: {
1128 double t = 2*y100 - 61;
1129 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;
1130 }
1131 case 31: {
1132 double t = 2*y100 - 63;
1133 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;
1134 }
1135 case 32: {
1136 double t = 2*y100 - 65;
1137 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;
1138 }
1139 case 33: {
1140 double t = 2*y100 - 67;
1141 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;
1142 }
1143 case 34: {
1144 double t = 2*y100 - 69;
1145 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;
1146 }
1147 case 35: {
1148 double t = 2*y100 - 71;
1149 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;
1150 }
1151 case 36: {
1152 double t = 2*y100 - 73;
1153 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;
1154 }
1155 case 37: {
1156 double t = 2*y100 - 75;
1157 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;
1158 }
1159 case 38: {
1160 double t = 2*y100 - 77;
1161 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;
1162 }
1163 case 39: {
1164 double t = 2*y100 - 79;
1165 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;
1166 }
1167 case 40: {
1168 double t = 2*y100 - 81;
1169 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;
1170 }
1171 case 41: {
1172 double t = 2*y100 - 83;
1173 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;
1174 }
1175 case 42: {
1176 double t = 2*y100 - 85;
1177 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;
1178 }
1179 case 43: {
1180 double t = 2*y100 - 87;
1181 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;
1182 }
1183 case 44: {
1184 double t = 2*y100 - 89;
1185 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;
1186 }
1187 case 45: {
1188 double t = 2*y100 - 91;
1189 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;
1190 }
1191 case 46: {
1192 double t = 2*y100 - 93;
1193 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;
1194 }
1195 case 47: {
1196 double t = 2*y100 - 95;
1197 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;
1198 }
1199 case 48: {
1200 double t = 2*y100 - 97;
1201 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;
1202 }
1203 case 49: {
1204 double t = 2*y100 - 99;
1205 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;
1206 }
1207 case 50: {
1208 double t = 2*y100 - 101;
1209 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;
1210 }
1211 case 51: {
1212 double t = 2*y100 - 103;
1213 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;
1214 }
1215 case 52: {
1216 double t = 2*y100 - 105;
1217 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;
1218 }
1219 case 53: {
1220 double t = 2*y100 - 107;
1221 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;
1222 }
1223 case 54: {
1224 double t = 2*y100 - 109;
1225 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;
1226 }
1227 case 55: {
1228 double t = 2*y100 - 111;
1229 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;
1230 }
1231 case 56: {
1232 double t = 2*y100 - 113;
1233 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;
1234 }
1235 case 57: {
1236 double t = 2*y100 - 115;
1237 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;
1238 }
1239 case 58: {
1240 double t = 2*y100 - 117;
1241 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;
1242 }
1243 case 59: {
1244 double t = 2*y100 - 119;
1245 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;
1246 }
1247 case 60: {
1248 double t = 2*y100 - 121;
1249 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;
1250 }
1251 case 61: {
1252 double t = 2*y100 - 123;
1253 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;
1254 }
1255 case 62: {
1256 double t = 2*y100 - 125;
1257 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;
1258 }
1259 case 63: {
1260 double t = 2*y100 - 127;
1261 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;
1262 }
1263 case 64: {
1264 double t = 2*y100 - 129;
1265 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;
1266 }
1267 case 65: {
1268 double t = 2*y100 - 131;
1269 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;
1270 }
1271 case 66: {
1272 double t = 2*y100 - 133;
1273 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;
1274 }
1275 case 67: {
1276 double t = 2*y100 - 135;
1277 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;
1278 }
1279 case 68: {
1280 double t = 2*y100 - 137;
1281 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;
1282 }
1283 case 69: {
1284 double t = 2*y100 - 139;
1285 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;
1286 }
1287 case 70: {
1288 double t = 2*y100 - 141;
1289 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;
1290 }
1291 case 71: {
1292 double t = 2*y100 - 143;
1293 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;
1294 }
1295 case 72: {
1296 double t = 2*y100 - 145;
1297 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;
1298 }
1299 case 73: {
1300 double t = 2*y100 - 147;
1301 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;
1302 }
1303 case 74: {
1304 double t = 2*y100 - 149;
1305 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;
1306 }
1307 case 75: {
1308 double t = 2*y100 - 151;
1309 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;
1310 }
1311 case 76: {
1312 double t = 2*y100 - 153;
1313 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;
1314 }
1315 case 77: {
1316 double t = 2*y100 - 155;
1317 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;
1318 }
1319 case 78: {
1320 double t = 2*y100 - 157;
1321 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;
1322 }
1323 case 79: {
1324 double t = 2*y100 - 159;
1325 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;
1326 }
1327 case 80: {
1328 double t = 2*y100 - 161;
1329 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;
1330 }
1331 case 81: {
1332 double t = 2*y100 - 163;
1333 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;
1334 }
1335 case 82: {
1336 double t = 2*y100 - 165;
1337 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;
1338 }
1339 case 83: {
1340 double t = 2*y100 - 167;
1341 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;
1342 }
1343 case 84: {
1344 double t = 2*y100 - 169;
1345 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;
1346 }
1347 case 85: {
1348 double t = 2*y100 - 171;
1349 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;
1350 }
1351 case 86: {
1352 double t = 2*y100 - 173;
1353 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;
1354 }
1355 case 87: {
1356 double t = 2*y100 - 175;
1357 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;
1358 }
1359 case 88: {
1360 double t = 2*y100 - 177;
1361 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;
1362 }
1363 case 89: {
1364 double t = 2*y100 - 179;
1365 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;
1366 }
1367 case 90: {
1368 double t = 2*y100 - 181;
1369 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;
1370 }
1371 case 91: {
1372 double t = 2*y100 - 183;
1373 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;
1374 }
1375 case 92: {
1376 double t = 2*y100 - 185;
1377 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;
1378 }
1379 case 93: {
1380 double t = 2*y100 - 187;
1381 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;
1382 }
1383 case 94: {
1384 double t = 2*y100 - 189;
1385 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;
1386 }
1387 case 95: {
1388 double t = 2*y100 - 191;
1389 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;
1390 }
1391 case 96: {
1392 double t = 2*y100 - 193;
1393 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;
1394 }
1395 case 97: {
1396 double t = 2*y100 - 195;
1397 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;
1398 }
1399 case 98: {
1400 double t = 2*y100 - 197;
1401 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;
1402 }
1403 case 99: {
1404 double t = 2*y100 - 199;
1405 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;
1406 }
1407  }
1408  // we only get here if y = 1, i.e. |x| < 4*eps, in which case
1409  // erfcx is within 1e-15 of 1..
1410  return 1.0;
1411 }
1412 
1413 double FADDEEVA_RE(erfcx)(double x)
1414 {
1415  if (x >= 0) {
1416  if (x > 50) { // continued-fraction expansion is faster
1417  const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1418  if (x > 5e7) // 1-term expansion, important to avoid overflow
1419  return ispi / x;
1420  /* 5-term expansion (rely on compiler for CSE), simplified from:
1421  ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */
1422  return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
1423  }
1424  return erfcx_y100(400/(4+x));
1425  }
1426  else
1427  return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x)
1428  : 2*exp(x*x) - erfcx_y100(400/(4-x)));
1429 }
1430 
1431 /////////////////////////////////////////////////////////////////////////
1432 /* Compute a scaled Dawson integral
1433  FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
1434  equivalent to the imaginary part w(x) for real x.
1435 
1436  Uses methods similar to the erfcx calculation above: continued fractions
1437  for large |x|, a lookup table of Chebyshev polynomials for smaller |x|,
1438  and finally a Taylor expansion for |x|<0.01.
1439 
1440  Steven G. Johnson, October 2012. */
1441 
1442 /* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).
1443 
1444  Uses a look-up table of 100 different Chebyshev polynomials
1445  for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
1446  with the help of Maple and a little shell script. This allows
1447  the Chebyshev polynomials to be of significantly lower degree (about 1/30)
1448  compared to fitting the whole [0,1] interval with a single polynomial. */
1449 static double w_im_y100(double y100, double x) {
1450  switch (static_cast<int> (y100)) {
1451  case 0: {
1452  double t = 2*y100 - 1;
1453  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;
1454  }
1455  case 1: {
1456  double t = 2*y100 - 3;
1457  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;
1458  }
1459  case 2: {
1460  double t = 2*y100 - 5;
1461  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;
1462  }
1463  case 3: {
1464  double t = 2*y100 - 7;
1465  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;
1466  }
1467  case 4: {
1468  double t = 2*y100 - 9;
1469  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;
1470  }
1471  case 5: {
1472  double t = 2*y100 - 11;
1473  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;
1474  }
1475  case 6: {
1476  double t = 2*y100 - 13;
1477  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;
1478  }
1479  case 7: {
1480  double t = 2*y100 - 15;
1481  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;
1482  }
1483  case 8: {
1484  double t = 2*y100 - 17;
1485  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;
1486  }
1487  case 9: {
1488  double t = 2*y100 - 19;
1489  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;
1490  }
1491  case 10: {
1492  double t = 2*y100 - 21;
1493  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;
1494  }
1495  case 11: {
1496  double t = 2*y100 - 23;
1497  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;
1498  }
1499  case 12: {
1500  double t = 2*y100 - 25;
1501  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;
1502  }
1503  case 13: {
1504  double t = 2*y100 - 27;
1505  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;
1506  }
1507  case 14: {
1508  double t = 2*y100 - 29;
1509  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;
1510  }
1511  case 15: {
1512  double t = 2*y100 - 31;
1513  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;
1514  }
1515  case 16: {
1516  double t = 2*y100 - 33;
1517  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;
1518  }
1519  case 17: {
1520  double t = 2*y100 - 35;
1521  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;
1522  }
1523  case 18: {
1524  double t = 2*y100 - 37;
1525  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;
1526  }
1527  case 19: {
1528  double t = 2*y100 - 39;
1529  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;
1530  }
1531  case 20: {
1532  double t = 2*y100 - 41;
1533  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;
1534  }
1535  case 21: {
1536  double t = 2*y100 - 43;
1537  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;
1538  }
1539  case 22: {
1540  double t = 2*y100 - 45;
1541  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;
1542  }
1543  case 23: {
1544  double t = 2*y100 - 47;
1545  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;
1546  }
1547  case 24: {
1548  double t = 2*y100 - 49;
1549  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;
1550  }
1551  case 25: {
1552  double t = 2*y100 - 51;
1553  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;
1554  }
1555  case 26: {
1556  double t = 2*y100 - 53;
1557  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;
1558  }
1559  case 27: {
1560  double t = 2*y100 - 55;
1561  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;
1562  }
1563  case 28: {
1564  double t = 2*y100 - 57;
1565  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;
1566  }
1567  case 29: {
1568  double t = 2*y100 - 59;
1569  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;
1570  }
1571  case 30: {
1572  double t = 2*y100 - 61;
1573  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;
1574  }
1575  case 31: {
1576  double t = 2*y100 - 63;
1577  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;
1578  }
1579  case 32: {
1580  double t = 2*y100 - 65;
1581  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;
1582  }
1583  case 33: {
1584  double t = 2*y100 - 67;
1585  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;
1586  }
1587  case 34: {
1588  double t = 2*y100 - 69;
1589  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;
1590  }
1591  case 35: {
1592  double t = 2*y100 - 71;
1593  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;
1594  }
1595  case 36: {
1596  double t = 2*y100 - 73;
1597  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;
1598  }
1599  case 37: {
1600  double t = 2*y100 - 75;
1601  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;
1602  }
1603  case 38: {
1604  double t = 2*y100 - 77;
1605  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;
1606  }
1607  case 39: {
1608  double t = 2*y100 - 79;
1609  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;
1610  }
1611  case 40: {
1612  double t = 2*y100 - 81;
1613  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;
1614  }
1615  case 41: {
1616  double t = 2*y100 - 83;
1617  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;
1618  }
1619  case 42: {
1620  double t = 2*y100 - 85;
1621  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;
1622  }
1623  case 43: {
1624  double t = 2*y100 - 87;
1625  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;
1626  }
1627  case 44: {
1628  double t = 2*y100 - 89;
1629  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;
1630  }
1631  case 45: {
1632  double t = 2*y100 - 91;
1633  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;
1634  }
1635  case 46: {
1636  double t = 2*y100 - 93;
1637  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;
1638  }
1639  case 47: {
1640  double t = 2*y100 - 95;
1641  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;
1642  }
1643  case 48: {
1644  double t = 2*y100 - 97;
1645  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;
1646  }
1647  case 49: {
1648  double t = 2*y100 - 99;
1649  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;
1650  }
1651  case 50: {
1652  double t = 2*y100 - 101;
1653  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;
1654  }
1655  case 51: {
1656  double t = 2*y100 - 103;
1657  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;
1658  }
1659  case 52: {
1660  double t = 2*y100 - 105;
1661  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;
1662  }
1663  case 53: {
1664  double t = 2*y100 - 107;
1665  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;
1666  }
1667  case 54: {
1668  double t = 2*y100 - 109;
1669  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;
1670  }
1671  case 55: {
1672  double t = 2*y100 - 111;
1673  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;
1674  }
1675  case 56: {
1676  double t = 2*y100 - 113;
1677  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;
1678  }
1679  case 57: {
1680  double t = 2*y100 - 115;
1681  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;
1682  }
1683  case 58: {
1684  double t = 2*y100 - 117;
1685  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;
1686  }
1687  case 59: {
1688  double t = 2*y100 - 119;
1689  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;
1690  }
1691  case 60: {
1692  double t = 2*y100 - 121;
1693  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;
1694  }
1695  case 61: {
1696  double t = 2*y100 - 123;
1697  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;
1698  }
1699  case 62: {
1700  double t = 2*y100 - 125;
1701  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;
1702  }
1703  case 63: {
1704  double t = 2*y100 - 127;
1705  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;
1706  }
1707  case 64: {
1708  double t = 2*y100 - 129;
1709  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;
1710  }
1711  case 65: {
1712  double t = 2*y100 - 131;
1713  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;
1714  }
1715  case 66: {
1716  double t = 2*y100 - 133;
1717  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;
1718  }
1719  case 67: {
1720  double t = 2*y100 - 135;
1721  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;
1722  }
1723  case 68: {
1724  double t = 2*y100 - 137;
1725  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;
1726  }
1727  case 69: {
1728  double t = 2*y100 - 139;
1729  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;
1730  }
1731  case 70: {
1732  double t = 2*y100 - 141;
1733  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;
1734  }
1735  case 71: {
1736  double t = 2*y100 - 143;
1737  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;
1738  }
1739  case 72: {
1740  double t = 2*y100 - 145;
1741  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;
1742  }
1743  case 73: {
1744  double t = 2*y100 - 147;
1745  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;
1746  }
1747  case 74: {
1748  double t = 2*y100 - 149;
1749  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;
1750  }
1751  case 75: {
1752  double t = 2*y100 - 151;
1753  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;
1754  }
1755  case 76: {
1756  double t = 2*y100 - 153;
1757  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;
1758  }
1759  case 77: {
1760  double t = 2*y100 - 155;
1761  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;
1762  }
1763  case 78: {
1764  double t = 2*y100 - 157;
1765  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;
1766  }
1767  case 79: {
1768  double t = 2*y100 - 159;
1769  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;
1770  }
1771  case 80: {
1772  double t = 2*y100 - 161;
1773  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;
1774  }
1775  case 81: {
1776  double t = 2*y100 - 163;
1777  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;
1778  }
1779  case 82: {
1780  double t = 2*y100 - 165;
1781  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;
1782  }
1783  case 83: {
1784  double t = 2*y100 - 167;
1785  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;
1786  }
1787  case 84: {
1788  double t = 2*y100 - 169;
1789  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;
1790  }
1791  case 85: {
1792  double t = 2*y100 - 171;
1793  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;
1794  }
1795  case 86: {
1796  double t = 2*y100 - 173;
1797  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;
1798  }
1799  case 87: {
1800  double t = 2*y100 - 175;
1801  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;
1802  }
1803  case 88: {
1804  double t = 2*y100 - 177;
1805  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;
1806  }
1807  case 89: {
1808  double t = 2*y100 - 179;
1809  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;
1810  }
1811  case 90: {
1812  double t = 2*y100 - 181;
1813  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;
1814  }
1815  case 91: {
1816  double t = 2*y100 - 183;
1817  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;
1818  }
1819  case 92: {
1820  double t = 2*y100 - 185;
1821  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;
1822  }
1823  case 93: {
1824  double t = 2*y100 - 187;
1825  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;
1826  }
1827  case 94: {
1828  double t = 2*y100 - 189;
1829  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;
1830  }
1831  case 95: {
1832  double t = 2*y100 - 191;
1833  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;
1834  }
1835  case 96: {
1836  double t = 2*y100 - 193;
1837  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;
1838  }
1839  case 97: case 98:
1840  case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
1841  // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9)
1842  double x2 = x*x;
1843  return x * (1.1283791670955125739
1844  - x2 * (0.75225277806367504925
1845  - x2 * (0.30090111122547001970
1846  - x2 * (0.085971746064420005629
1847  - x2 * 0.016931216931216931217))));
1848  }
1849  }
1850  /* Since 0 <= y100 < 101, this is only reached if x is NaN,
1851  in which case we should return NaN. */
1852  return NaN;
1853 }
1854 
1855 double FADDEEVA(w_im)(double x)
1856 {
1857  if (x >= 0) {
1858  if (x > 45) { // continued-fraction expansion is faster
1859  const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1860  if (x > 5e7) // 1-term expansion, important to avoid overflow
1861  return ispi / x;
1862  /* 5-term expansion (rely on compiler for CSE), simplified from:
1863  ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1864  return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1865  }
1866  return w_im_y100(100/(1+x), x);
1867  }
1868  else { // = -FADDEEVA(w_im)(-x)
1869  if (x < -45) { // continued-fraction expansion is faster
1870  const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1871  if (x < -5e7) // 1-term expansion, important to avoid overflow
1872  return ispi / x;
1873  /* 5-term expansion (rely on compiler for CSE), simplified from:
1874  ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1875  return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1876  }
1877  return -w_im_y100(100/(1-x), -x);
1878  }
1879 }
1880 
1881 /////////////////////////////////////////////////////////////////////////
1882 
1883 // Compile with -DTEST_FADDEEVA to compile a little test program
1884 #if defined (TEST_FADDEEVA)
1885 
1886 #if defined (__cplusplus)
1887 # include <cstdio>
1888 #else
1889 # include <stdio.h>
1890 #endif
1891 
1892 // compute relative error |b-a|/|a|, handling case of NaN and Inf,
1893 static double relerr(double a, double b) {
1894  if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) {
1895  if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) ||
1896  (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) ||
1897  (isinf(a) && isinf(b) && a*b < 0))
1898  return Inf; // "infinite" error
1899  return 0; // matching infinity/nan results counted as zero error
1900  }
1901  if (a == 0)
1902  return b == 0 ? 0 : Inf;
1903  else
1904  return fabs((b-a) / a);
1905 }
1906 
1907 int main(void) {
1908  double errmax_all = 0;
1909  {
1910  printf("############# w(z) tests #############\n");
1911 #define NTST 57 // define instead of const for C compatibility
1912  cmplx z[NTST] = {
1913  C(624.2,-0.26123),
1914  C(-0.4,3.),
1915  C(0.6,2.),
1916  C(-1.,1.),
1917  C(-1.,-9.),
1918  C(-1.,9.),
1919  C(-0.0000000234545,1.1234),
1920  C(-3.,5.1),
1921  C(-53,30.1),
1922  C(0.0,0.12345),
1923  C(11,1),
1924  C(-22,-2),
1925  C(9,-28),
1926  C(21,-33),
1927  C(1e5,1e5),
1928  C(1e14,1e14),
1929  C(-3001,-1000),
1930  C(1e160,-1e159),
1931  C(-6.01,0.01),
1932  C(-0.7,-0.7),
1933  C(2.611780000000000e+01, 4.540909610972489e+03),
1934  C(0.8e7,0.3e7),
1935  C(-20,-19.8081),
1936  C(1e-16,-1.1e-16),
1937  C(2.3e-8,1.3e-8),
1938  C(6.3,-1e-13),
1939  C(6.3,1e-20),
1940  C(1e-20,6.3),
1941  C(1e-20,16.3),
1942  C(9,1e-300),
1943  C(6.01,0.11),
1944  C(8.01,1.01e-10),
1945  C(28.01,1e-300),
1946  C(10.01,1e-200),
1947  C(10.01,-1e-200),
1948  C(10.01,0.99e-10),
1949  C(10.01,-0.99e-10),
1950  C(1e-20,7.01),
1951  C(-1,7.01),
1952  C(5.99,7.01),
1953  C(1,0),
1954  C(55,0),
1955  C(-0.1,0),
1956  C(1e-20,0),
1957  C(0,5e-14),
1958  C(0,51),
1959  C(Inf,0),
1960  C(-Inf,0),
1961  C(0,Inf),
1962  C(0,-Inf),
1963  C(Inf,Inf),
1964  C(Inf,-Inf),
1965  C(NaN,NaN),
1966  C(NaN,0),
1967  C(0,NaN),
1968  C(NaN,Inf),
1969  C(Inf,NaN)
1970  };
1971  cmplx w[NTST] = { /* w(z), computed with WolframAlpha
1972  ... note that WolframAlpha is problematic
1973  some of the above inputs, so I had to
1974  use the continued-fraction expansion
1975  in WolframAlpha in some cases, or switch
1976  to Maple */
1977  C(-3.78270245518980507452677445620103199303131110e-7,
1978  0.000903861276433172057331093754199933411710053155),
1979  C(0.1764906227004816847297495349730234591778719532788,
1980  -0.02146550539468457616788719893991501311573031095617),
1981  C(0.2410250715772692146133539023007113781272362309451,
1982  0.06087579663428089745895459735240964093522265589350),
1983  C(0.30474420525691259245713884106959496013413834051768,
1984  -0.20821893820283162728743734725471561394145872072738),
1985  C(7.317131068972378096865595229600561710140617977e34,
1986  8.321873499714402777186848353320412813066170427e34),
1987  C(0.0615698507236323685519612934241429530190806818395,
1988  -0.00676005783716575013073036218018565206070072304635),
1989  C(0.3960793007699874918961319170187598400134746631,
1990  -5.593152259116644920546186222529802777409274656e-9),
1991  C(0.08217199226739447943295069917990417630675021771804,
1992  -0.04701291087643609891018366143118110965272615832184),
1993  C(0.00457246000350281640952328010227885008541748668738,
1994  -0.00804900791411691821818731763401840373998654987934),
1995  C(0.8746342859608052666092782112565360755791467973338452,
1996  0.),
1997  C(0.00468190164965444174367477874864366058339647648741,
1998  0.0510735563901306197993676329845149741675029197050),
1999  C(-0.0023193175200187620902125853834909543869428763219,
2000  -0.025460054739731556004902057663500272721780776336),
2001  C(9.11463368405637174660562096516414499772662584e304,
2002  3.97101807145263333769664875189354358563218932e305),
2003  C(-4.4927207857715598976165541011143706155432296e281,
2004  -2.8019591213423077494444700357168707775769028e281),
2005  C(2.820947917809305132678577516325951485807107151e-6,
2006  2.820947917668257736791638444590253942253354058e-6),
2007  C(2.82094791773878143474039725787438662716372268e-15,
2008  2.82094791773878143474039725773333923127678361e-15),
2009  C(-0.0000563851289696244350147899376081488003110150498,
2010  -0.000169211755126812174631861529808288295454992688),
2011  C(-5.586035480670854326218608431294778077663867e-162,
2012  5.586035480670854326218608431294778077663867e-161),
2013  C(0.00016318325137140451888255634399123461580248456,
2014  -0.095232456573009287370728788146686162555021209999),
2015  C(0.69504753678406939989115375989939096800793577783885,
2016  -1.8916411171103639136680830887017670616339912024317),
2017  C(0.0001242418269653279656612334210746733213167234822,
2018  7.145975826320186888508563111992099992116786763e-7),
2019  C(2.318587329648353318615800865959225429377529825e-8,
2020  6.182899545728857485721417893323317843200933380e-8),
2021  C(-0.0133426877243506022053521927604277115767311800303,
2022  -0.0148087097143220769493341484176979826888871576145),
2023  C(1.00000000000000012412170838050638522857747934,
2024  1.12837916709551279389615890312156495593616433e-16),
2025  C(0.9999999853310704677583504063775310832036830015,
2026  2.595272024519678881897196435157270184030360773e-8),
2027  C(-1.4731421795638279504242963027196663601154624e-15,
2028  0.090727659684127365236479098488823462473074709),
2029  C(5.79246077884410284575834156425396800754409308e-18,
2030  0.0907276596841273652364790985059772809093822374),
2031  C(0.0884658993528521953466533278764830881245144368,
2032  1.37088352495749125283269718778582613192166760e-22),
2033  C(0.0345480845419190424370085249304184266813447878,
2034  2.11161102895179044968099038990446187626075258e-23),
2035  C(6.63967719958073440070225527042829242391918213e-36,
2036  0.0630820900592582863713653132559743161572639353),
2037  C(0.00179435233208702644891092397579091030658500743634,
2038  0.0951983814805270647939647438459699953990788064762),
2039  C(9.09760377102097999924241322094863528771095448e-13,
2040  0.0709979210725138550986782242355007611074966717),
2041  C(7.2049510279742166460047102593255688682910274423e-304,
2042  0.0201552956479526953866611812593266285000876784321),
2043  C(3.04543604652250734193622967873276113872279682e-44,
2044  0.0566481651760675042930042117726713294607499165),
2045  C(3.04543604652250734193622967873276113872279682e-44,
2046  0.0566481651760675042930042117726713294607499165),
2047  C(0.5659928732065273429286988428080855057102069081e-12,
2048  0.056648165176067504292998527162143030538756683302),
2049  C(-0.56599287320652734292869884280802459698927645e-12,
2050  0.0566481651760675042929985271621430305387566833029),
2051  C(0.0796884251721652215687859778119964009569455462,
2052  1.11474461817561675017794941973556302717225126e-22),
2053  C(0.07817195821247357458545539935996687005781943386550,
2054  -0.01093913670103576690766705513142246633056714279654),
2055  C(0.04670032980990449912809326141164730850466208439937,
2056  0.03944038961933534137558064191650437353429669886545),
2057  C(0.36787944117144232159552377016146086744581113103176,
2058  0.60715770584139372911503823580074492116122092866515),
2059  C(0,
2060  0.010259688805536830986089913987516716056946786526145),
2061  C(0.99004983374916805357390597718003655777207908125383,
2062  -0.11208866436449538036721343053869621153527769495574),
2063  C(0.99999999999999999999999999999999999999990000,
2064  1.12837916709551257389615890312154517168802603e-20),
2065  C(0.999999999999943581041645226871305192054749891144158,
2066  0),
2067  C(0.0110604154853277201542582159216317923453996211744250,
2068  0),
2069  C(0,0),
2070  C(0,0),
2071  C(0,0),
2072  C(Inf,0),
2073  C(0,0),
2074  C(NaN,NaN),
2075  C(NaN,NaN),
2076  C(NaN,NaN),
2077  C(NaN,0),
2078  C(NaN,NaN),
2079  C(NaN,NaN)
2080  };
2081  double errmax = 0;
2082  for (int i = 0; i < NTST; ++i) {
2083  cmplx fw = FADDEEVA(w)(z[i],0.);
2084  double re_err = relerr(creal(w[i]), creal(fw));
2085  double im_err = relerr(cimag(w[i]), cimag(fw));
2086  printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n",
2087  creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]),
2088  re_err, im_err);
2089  if (re_err > errmax) errmax = re_err;
2090  if (im_err > errmax) errmax = im_err;
2091  }
2092  if (errmax > 1e-13) {
2093  printf("FAILURE -- relative error %g too large!\n", errmax);
2094  return 1;
2095  }
2096  printf("SUCCESS (max relative error = %g)\n", errmax);
2097  if (errmax > errmax_all) errmax_all = errmax;
2098  }
2099  {
2100 #undef NTST
2101 #define NTST 41 // define instead of const for C compatibility
2102  cmplx z[NTST] = {
2103  C(1,2),
2104  C(-1,2),
2105  C(1,-2),
2106  C(-1,-2),
2107  C(9,-28),
2108  C(21,-33),
2109  C(1e3,1e3),
2110  C(-3001,-1000),
2111  C(1e160,-1e159),
2112  C(5.1e-3, 1e-8),
2113  C(-4.9e-3, 4.95e-3),
2114  C(4.9e-3, 0.5),
2115  C(4.9e-4, -0.5e1),
2116  C(-4.9e-5, -0.5e2),
2117  C(5.1e-3, 0.5),
2118  C(5.1e-4, -0.5e1),
2119  C(-5.1e-5, -0.5e2),
2120  C(1e-6,2e-6),
2121  C(0,2e-6),
2122  C(0,2),
2123  C(0,20),
2124  C(0,200),
2125  C(Inf,0),
2126  C(-Inf,0),
2127  C(0,Inf),
2128  C(0,-Inf),
2129  C(Inf,Inf),
2130  C(Inf,-Inf),
2131  C(NaN,NaN),
2132  C(NaN,0),
2133  C(0,NaN),
2134  C(NaN,Inf),
2135  C(Inf,NaN),
2136  C(1e-3,NaN),
2137  C(7e-2,7e-2),
2138  C(7e-2,-7e-4),
2139  C(-9e-2,7e-4),
2140  C(-9e-2,9e-2),
2141  C(-7e-4,9e-2),
2142  C(7e-2,0.9e-2),
2143  C(7e-2,1.1e-2)
2144  };
2145  cmplx w[NTST] = { // erf(z[i]), evaluated with Maple
2146  C(-0.5366435657785650339917955593141927494421,
2147  -5.049143703447034669543036958614140565553),
2148  C(0.5366435657785650339917955593141927494421,
2149  -5.049143703447034669543036958614140565553),
2150  C(-0.5366435657785650339917955593141927494421,
2151  5.049143703447034669543036958614140565553),
2152  C(0.5366435657785650339917955593141927494421,
2153  5.049143703447034669543036958614140565553),
2154  C(0.3359473673830576996788000505817956637777e304,
2155  -0.1999896139679880888755589794455069208455e304),
2156  C(0.3584459971462946066523939204836760283645e278,
2157  0.3818954885257184373734213077678011282505e280),
2158  C(0.9996020422657148639102150147542224526887,
2159  0.00002801044116908227889681753993542916894856),
2160  C(-1, 0),
2161  C(1, 0),
2162  C(0.005754683859034800134412990541076554934877,
2163  0.1128349818335058741511924929801267822634e-7),
2164  C(-0.005529149142341821193633460286828381876955,
2165  0.005585388387864706679609092447916333443570),
2166  C(0.007099365669981359632319829148438283865814,
2167  0.6149347012854211635026981277569074001219),
2168  C(0.3981176338702323417718189922039863062440e8,
2169  -0.8298176341665249121085423917575122140650e10),
2170  C(-Inf,
2171  -Inf),
2172  C(0.007389128308257135427153919483147229573895,
2173  0.6149332524601658796226417164791221815139),
2174  C(0.4143671923267934479245651547534414976991e8,
2175  -0.8298168216818314211557046346850921446950e10),
2176  C(-Inf,
2177  -Inf),
2178  C(0.1128379167099649964175513742247082845155e-5,
2179  0.2256758334191777400570377193451519478895e-5),
2180  C(0,
2181  0.2256758334194034158904576117253481476197e-5),
2182  C(0,
2183  18.56480241457555259870429191324101719886),
2184  C(0,
2185  0.1474797539628786202447733153131835124599e173),
2186  C(0,
2187  Inf),
2188  C(1,0),
2189  C(-1,0),
2190  C(0,Inf),
2191  C(0,-Inf),
2192  C(NaN,NaN),
2193  C(NaN,NaN),
2194  C(NaN,NaN),
2195  C(NaN,0),
2196  C(0,NaN),
2197  C(NaN,NaN),
2198  C(NaN,NaN),
2199  C(NaN,NaN),
2200  C(0.07924380404615782687930591956705225541145,
2201  0.07872776218046681145537914954027729115247),
2202  C(0.07885775828512276968931773651224684454495,
2203  -0.0007860046704118224342390725280161272277506),
2204  C(-0.1012806432747198859687963080684978759881,
2205  0.0007834934747022035607566216654982820299469),
2206  C(-0.1020998418798097910247132140051062512527,
2207  0.1010030778892310851309082083238896270340),
2208  C(-0.0007962891763147907785684591823889484764272,
2209  0.1018289385936278171741809237435404896152),
2210  C(0.07886408666470478681566329888615410479530,
2211  0.01010604288780868961492224347707949372245),
2212  C(0.07886723099940260286824654364807981336591,
2213  0.01235199327873258197931147306290916629654)
2214  };
2215 #define TST(f,isc) \
2216  printf("############# " #f "(z) tests #############\n"); \
2217  double errmax = 0; \
2218  for (int i = 0; i < NTST; ++i) { \
2219  cmplx fw = FADDEEVA(f)(z[i],0.); \
2220  double re_err = relerr(creal(w[i]), creal(fw)); \
2221  double im_err = relerr(cimag(w[i]), cimag(fw)); \
2222  printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \
2223  creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \
2224  re_err, im_err); \
2225  if (re_err > errmax) errmax = re_err; \
2226  if (im_err > errmax) errmax = im_err; \
2227  } \
2228  if (errmax > 1e-13) { \
2229  printf("FAILURE -- relative error %g too large!\n", errmax); \
2230  return 1; \
2231  } \
2232  printf("Checking " #f "(x) special case...\n"); \
2233  for (int i = 0; i < 10000; ++i) { \
2234  double x = pow(10., -300. + i * 600. / (10000 - 1)); \
2235  double re_err = relerr(FADDEEVA_RE(f)(x), \
2236  creal(FADDEEVA(f)(C(x,x*isc),0.))); \
2237  if (re_err > errmax) errmax = re_err; \
2238  re_err = relerr(FADDEEVA_RE(f)(-x), \
2239  creal(FADDEEVA(f)(C(-x,x*isc),0.))); \
2240  if (re_err > errmax) errmax = re_err; \
2241  } \
2242  { \
2243  double re_err = relerr(FADDEEVA_RE(f)(Inf), \
2244  creal(FADDEEVA(f)(C(Inf,0.),0.))); \
2245  if (re_err > errmax) errmax = re_err; \
2246  re_err = relerr(FADDEEVA_RE(f)(-Inf), \
2247  creal(FADDEEVA(f)(C(-Inf,0.),0.))); \
2248  if (re_err > errmax) errmax = re_err; \
2249  re_err = relerr(FADDEEVA_RE(f)(NaN), \
2250  creal(FADDEEVA(f)(C(NaN,0.),0.))); \
2251  if (re_err > errmax) errmax = re_err; \
2252  } \
2253  if (errmax > 1e-13) { \
2254  printf("FAILURE -- relative error %g too large!\n", errmax); \
2255  return 1; \
2256  } \
2257  printf("SUCCESS (max relative error = %g)\n", errmax); \
2258  if (errmax > errmax_all) errmax_all = errmax
2259 
2260  TST(erf, 1e-20);
2261  }
2262  {
2263  // since erfi just calls through to erf, just one test should
2264  // be sufficient to make sure I didn't screw up the signs or something
2265 #undef NTST
2266 #define NTST 1 // define instead of const for C compatibility
2267  cmplx z[NTST] = { C(1.234,0.5678) };
2268  cmplx w[NTST] = { // erfi(z[i]), computed with Maple
2269  C(1.081032284405373149432716643834106923212,
2270  1.926775520840916645838949402886591180834)
2271  };
2272  TST(erfi, 0);
2273  }
2274  {
2275  // since erfcx just calls through to w, just one test should
2276  // be sufficient to make sure I didn't screw up the signs or something
2277 #undef NTST
2278 #define NTST 1 // define instead of const for C compatibility
2279  cmplx z[NTST] = { C(1.234,0.5678) };
2280  cmplx w[NTST] = { // erfcx(z[i]), computed with Maple
2281  C(0.3382187479799972294747793561190487832579,
2282  -0.1116077470811648467464927471872945833154)
2283  };
2284  TST(erfcx, 0);
2285  }
2286  {
2287 #undef NTST
2288 #define NTST 30 // define instead of const for C compatibility
2289  cmplx z[NTST] = {
2290  C(1,2),
2291  C(-1,2),
2292  C(1,-2),
2293  C(-1,-2),
2294  C(9,-28),
2295  C(21,-33),
2296  C(1e3,1e3),
2297  C(-3001,-1000),
2298  C(1e160,-1e159),
2299  C(5.1e-3, 1e-8),
2300  C(0,2e-6),
2301  C(0,2),
2302  C(0,20),
2303  C(0,200),
2304  C(2e-6,0),
2305  C(2,0),
2306  C(20,0),
2307  C(200,0),
2308  C(Inf,0),
2309  C(-Inf,0),
2310  C(0,Inf),
2311  C(0,-Inf),
2312  C(Inf,Inf),
2313  C(Inf,-Inf),
2314  C(NaN,NaN),
2315  C(NaN,0),
2316  C(0,NaN),
2317  C(NaN,Inf),
2318  C(Inf,NaN),
2319  C(88,0)
2320  };
2321  cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple
2322  C(1.536643565778565033991795559314192749442,
2323  5.049143703447034669543036958614140565553),
2324  C(0.4633564342214349660082044406858072505579,
2325  5.049143703447034669543036958614140565553),
2326  C(1.536643565778565033991795559314192749442,
2327  -5.049143703447034669543036958614140565553),
2328  C(0.4633564342214349660082044406858072505579,
2329  -5.049143703447034669543036958614140565553),
2330  C(-0.3359473673830576996788000505817956637777e304,
2331  0.1999896139679880888755589794455069208455e304),
2332  C(-0.3584459971462946066523939204836760283645e278,
2333  -0.3818954885257184373734213077678011282505e280),
2334  C(0.0003979577342851360897849852457775473112748,
2335  -0.00002801044116908227889681753993542916894856),
2336  C(2, 0),
2337  C(0, 0),
2338  C(0.9942453161409651998655870094589234450651,
2339  -0.1128349818335058741511924929801267822634e-7),
2340  C(1,
2341  -0.2256758334194034158904576117253481476197e-5),
2342  C(1,
2343  -18.56480241457555259870429191324101719886),
2344  C(1,
2345  -0.1474797539628786202447733153131835124599e173),
2346  C(1, -Inf),
2347  C(0.9999977432416658119838633199332831406314,
2348  0),
2349  C(0.004677734981047265837930743632747071389108,
2350  0),
2351  C(0.5395865611607900928934999167905345604088e-175,
2352  0),
2353  C(0, 0),
2354  C(0, 0),
2355  C(2, 0),
2356  C(1, -Inf),
2357  C(1, Inf),
2358  C(NaN, NaN),
2359  C(NaN, NaN),
2360  C(NaN, NaN),
2361  C(NaN, 0),
2362  C(1, NaN),
2363  C(NaN, NaN),
2364  C(NaN, NaN),
2365  C(0,0)
2366  };
2367  TST(erfc, 1e-20);
2368  }
2369  {
2370 #undef NTST
2371 #define NTST 48 // define instead of const for C compatibility
2372  cmplx z[NTST] = {
2373  C(2,1),
2374  C(-2,1),
2375  C(2,-1),
2376  C(-2,-1),
2377  C(-28,9),
2378  C(33,-21),
2379  C(1e3,1e3),
2380  C(-1000,-3001),
2381  C(1e-8, 5.1e-3),
2382  C(4.95e-3, -4.9e-3),
2383  C(5.1e-3, 5.1e-3),
2384  C(0.5, 4.9e-3),
2385  C(-0.5e1, 4.9e-4),
2386  C(-0.5e2, -4.9e-5),
2387  C(0.5e3, 4.9e-6),
2388  C(0.5, 5.1e-3),
2389  C(-0.5e1, 5.1e-4),
2390  C(-0.5e2, -5.1e-5),
2391  C(1e-6,2e-6),
2392  C(2e-6,0),
2393  C(2,0),
2394  C(20,0),
2395  C(200,0),
2396  C(0,4.9e-3),
2397  C(0,-5.1e-3),
2398  C(0,2e-6),
2399  C(0,-2),
2400  C(0,20),
2401  C(0,-200),
2402  C(Inf,0),
2403  C(-Inf,0),
2404  C(0,Inf),
2405  C(0,-Inf),
2406  C(Inf,Inf),
2407  C(Inf,-Inf),
2408  C(NaN,NaN),
2409  C(NaN,0),
2410  C(0,NaN),
2411  C(NaN,Inf),
2412  C(Inf,NaN),
2413  C(39, 6.4e-5),
2414  C(41, 6.09e-5),
2415  C(4.9e7, 5e-11),
2416  C(5.1e7, 4.8e-11),
2417  C(1e9, 2.4e-12),
2418  C(1e11, 2.4e-14),
2419  C(1e13, 2.4e-16),
2420  C(1e300, 2.4e-303)
2421  };
2422  cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple
2423  C(0.1635394094345355614904345232875688576839,
2424  -0.1531245755371229803585918112683241066853),
2425  C(-0.1635394094345355614904345232875688576839,
2426  -0.1531245755371229803585918112683241066853),
2427  C(0.1635394094345355614904345232875688576839,
2428  0.1531245755371229803585918112683241066853),
2429  C(-0.1635394094345355614904345232875688576839,
2430  0.1531245755371229803585918112683241066853),
2431  C(-0.01619082256681596362895875232699626384420,
2432  -0.005210224203359059109181555401330902819419),
2433  C(0.01078377080978103125464543240346760257008,
2434  0.006866888783433775382193630944275682670599),
2435  C(-0.5808616819196736225612296471081337245459,
2436  0.6688593905505562263387760667171706325749),
2437  C(Inf,
2438  -Inf),
2439  C(0.1000052020902036118082966385855563526705e-7,
2440  0.005100088434920073153418834680320146441685),
2441  C(0.004950156837581592745389973960217444687524,
2442  -0.004899838305155226382584756154100963570500),
2443  C(0.005100176864319675957314822982399286703798,
2444  0.005099823128319785355949825238269336481254),
2445  C(0.4244534840871830045021143490355372016428,
2446  0.002820278933186814021399602648373095266538),
2447  C(-0.1021340733271046543881236523269967674156,
2448  -0.00001045696456072005761498961861088944159916),
2449  C(-0.01000200120119206748855061636187197886859,
2450  0.9805885888237419500266621041508714123763e-8),
2451  C(0.001000002000012000023960527532953151819595,
2452  -0.9800058800588007290937355024646722133204e-11),
2453  C(0.4244549085628511778373438768121222815752,
2454  0.002935393851311701428647152230552122898291),
2455  C(-0.1021340732357117208743299813648493928105,
2456  -0.00001088377943049851799938998805451564893540),
2457  C(-0.01000200120119126652710792390331206563616,
2458  0.1020612612857282306892368985525393707486e-7),
2459  C(0.1000000000007333333333344266666666664457e-5,
2460  0.2000000000001333333333323199999999978819e-5),
2461  C(0.1999999999994666666666675199999999990248e-5,
2462  0),
2463  C(0.3013403889237919660346644392864226952119,
2464  0),
2465  C(0.02503136792640367194699495234782353186858,
2466  0),
2467  C(0.002500031251171948248596912483183760683918,
2468  0),
2469  C(0,0.004900078433419939164774792850907128053308),
2470  C(0,-0.005100088434920074173454208832365950009419),
2471  C(0,0.2000000000005333333333341866666666676419e-5),
2472  C(0,-48.16001211429122974789822893525016528191),
2473  C(0,0.4627407029504443513654142715903005954668e174),
2474  C(0,-Inf),
2475  C(0,0),
2476  C(-0,0),
2477  C(0, Inf),
2478  C(0, -Inf),
2479  C(NaN, NaN),
2480  C(NaN, NaN),
2481  C(NaN, NaN),
2482  C(NaN, 0),
2483  C(0, NaN),
2484  C(NaN, NaN),
2485  C(NaN, NaN),
2486  C(0.01282473148489433743567240624939698290584,
2487  -0.2105957276516618621447832572909153498104e-7),
2488  C(0.01219875253423634378984109995893708152885,
2489  -0.1813040560401824664088425926165834355953e-7),
2490  C(0.1020408163265306334945473399689037886997e-7,
2491  -0.1041232819658476285651490827866174985330e-25),
2492  C(0.9803921568627452865036825956835185367356e-8,
2493  -0.9227220299884665067601095648451913375754e-26),
2494  C(0.5000000000000000002500000000000000003750e-9,
2495  -0.1200000000000000001800000188712838420241e-29),
2496  C(5.00000000000000000000025000000000000000000003e-12,
2497  -1.20000000000000000000018000000000000000000004e-36),
2498  C(5.00000000000000000000000002500000000000000000e-14,
2499  -1.20000000000000000000000001800000000000000000e-42),
2500  C(5e-301, 0)
2501  };
2502  TST(Dawson, 1e-20);
2503  }
2504  printf("#####################################\n");
2505  printf("SUCCESS (max relative error = %g)\n", errmax_all);
2506 }
2507 
2508 #endif
int main(int argc, char **argv)
Definition: main-cli.cc:82
#define C(a, b)
Definition: Faddeeva.cc:246
function erf(X)
Definition: erf.f:2
nd group nd example oindent or xample printf("%s\n", nthargout(2,"system","cmd"))
std::complex< double > erfi(std::complex< double > z, double relerr=0)
bool isnan(double x)
Definition: lo-mappers.cc:347
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
std::complex< double > Dawson(std::complex< double > z, double relerr=0)
static double sinc(double x, double sinx)
Definition: Faddeeva.cc:600
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
Definition: ov-usr-fcn.cc:935
i e
Definition: data.cc:2724
#define NaN
Definition: Faddeeva.cc:248
const octave_base_value & a2
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
static double w_im_y100(double y100, double x)
Definition: Faddeeva.cc:1449
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
Definition: utils.cc:941
double copysign(double x, double y)
Definition: lo-mappers.cc:318
std::complex< double > w(std::complex< double > z, double relerr=0)
#define FADDEEVA(name)
Definition: Faddeeva.cc:219
double complex cmplx
Definition: Faddeeva.cc:217
std::complex< double > erfcx(std::complex< double > z, double relerr=0)
static double sinh_taylor(double x)
Definition: Faddeeva.cc:605
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
#define FADDEEVA_RE(name)
Definition: Faddeeva.cc:220
static double sqr(double x)
Definition: Faddeeva.cc:610
bool isinf(double x)
Definition: lo-mappers.cc:387
double w_im(double x)
static cmplx cpolar(double r, double t)
Definition: Faddeeva.cc:251
#define Inf
Definition: Faddeeva.cc:247
static double wi[256]
Definition: randmtzig.cc:440
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
the element is set to zero In other the statement xample y
Definition: data.cc:5342
b
Definition: cellfun.cc:398
double floor(double x)
Definition: lo-mappers.cc:330
static const double expa2n2[]
Definition: Faddeeva.cc:614
static double erfcx_y100(double y100)
Definition: Faddeeva.cc:1004
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
static const double pi
Definition: lo-specfun.cc:3610
std::complex< double > erfc(std::complex< double > z, double relerr=0)