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 12 May 2015: Bugfix for systems lacking copysign function.
144 */
145
146 /////////////////////////////////////////////////////////////////////////
147 /* If this file is compiled as a part of a larger project,
148 support using an autoconf-style config.h header file
149 (with various "HAVE_*" #defines to indicate features)
150 if HAVE_CONFIG_H is #defined (in GNU autotools style). */
151
152 #ifdef HAVE_CONFIG_H
153 # include "libmesh_config.h"
154 #endif
155
156 /////////////////////////////////////////////////////////////////////////
157 // macros to allow us to use either C++ or C (with C99 features)
158
159 #ifdef __cplusplus
160
161 # include "Faddeeva.hh"
162
163 # include <cfloat>
164 # include <cmath>
165 # include <limits>
166 using namespace std;
167
168 // use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS)
169 # define Inf numeric_limits<double>::infinity()
170 # define NaN numeric_limits<double>::quiet_NaN()
171
172 typedef 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 (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF))
my_isnan(double x)187 static inline bool my_isnan(double x) { return x != x; }
188 # define isnan my_isnan
my_isinf(double x)189 static inline bool my_isinf(double x) { return 1/x == 0.; }
190 # define isinf my_isinf
191 # elif (__cplusplus >= 201103L)
192 // g++ gets confused between the C and C++ isnan/isinf functions
193 # define isnan std::isnan
194 # define isinf std::isinf
195 # endif
196
197 // copysign was introduced in C++11 (and is also in POSIX and C99)
198 # if defined(_WIN32) || defined(__WIN32__)
199 # define copysign _copysign // of course MS had to be different
200 # elif defined(GNULIB_NAMESPACE) // we are using using gnulib <cmath>
201 # define copysign GNULIB_NAMESPACE::copysign
202 # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX)
my_copysign(double x,double y)203 static inline double my_copysign(double x, double y) { return x<0 != y<0 ? -x : x; }
204 # define copysign my_copysign
205 # endif
206
207 // If we are using the gnulib <cmath> (e.g. in the GNU Octave sources),
208 // gnulib generates a link warning if we use ::floor instead of gnulib::floor.
209 // This warning is completely innocuous because the only difference between
210 // gnulib::floor and the system ::floor (and only on ancient OSF systems)
211 // has to do with floor(-0), which doesn't occur in the usage below, but
212 // the Octave developers prefer that we silence the warning.
213 # ifdef GNULIB_NAMESPACE
214 # define floor GNULIB_NAMESPACE::floor
215 # endif
216
217 #else // !__cplusplus, i.e. pure C (requires C99 features)
218
219 # include "Faddeeva.h"
220
221 # define _GNU_SOURCE // enable GNU libc NAN extension if possible
222
223 # include <float.h>
224 # include <math.h>
225
226 typedef double complex cmplx;
227
228 # define FADDEEVA(name) Faddeeva_ ## name
229 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
230
231 /* Constructing complex numbers like 0+i*NaN is problematic in C99
232 without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
233 I is a complex (rather than imaginary) constant. For some reason,
234 however, it works fine in (pre-4.7) gcc if I define Inf and NaN as
235 1/0 and 0/0 (and only if I compile with optimization -O1 or more),
236 but not if I use the INFINITY or NAN macros. */
237
238 /* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro
239 may not be defined unless we are using a recent (2012) version of
240 glibc and compile with -std=c11... note that icc lies about being
241 gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
242 # if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER))
243 # define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
244 # endif
245
246 # ifdef CMPLX // C11
247 # define C(a,b) CMPLX(a,b)
248 # define Inf INFINITY // C99 infinity
249 # ifdef NAN // GNU libc extension
250 # define NaN NAN
251 # else
252 # define NaN (0./0.) // NaN
253 # endif
254 # else
255 # define C(a,b) ((a) + I*(b))
256 # define Inf (1./0.)
257 # define NaN (0./0.)
258 # endif
259
cpolar(double r,double t)260 static inline cmplx cpolar(double r, double t)
261 {
262 if (r == 0.0 && !isnan(t))
263 return 0.0;
264 else
265 return C(r * cos(t), r * sin(t));
266 }
267
268 #endif // !__cplusplus, i.e. pure C (requires C99 features)
269
270 /////////////////////////////////////////////////////////////////////////
271 // Auxiliary routines to compute other special functions based on w(z)
272
273 // compute erfcx(z) = exp(z^2) erfz(z)
FADDEEVA(erfcx)274 cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
275 {
276 return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
277 }
278
279 // compute the error function erf(x)
FADDEEVA_RE(erf)280 double FADDEEVA_RE(erf)(double x)
281 {
282 #if !defined(__cplusplus)
283 return erf(x); // C99 supplies erf in math.h
284 #elif (__cplusplus >= 201103L) || defined(HAVE_ERF)
285 return ::erf(x); // C++11 supplies std::erf in cmath
286 #else
287 double mx2 = -x*x;
288 if (mx2 < -750) // underflow
289 return (x >= 0 ? 1.0 : -1.0);
290
291 if (x >= 0) {
292 if (x < 8e-2) goto taylor;
293 return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
294 }
295 else { // x < 0
296 if (x > -8e-2) goto taylor;
297 return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0;
298 }
299
300 // Use Taylor series for small |x|, to avoid cancellation inaccuracy
301 // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...)
302 taylor:
303 return x * (1.1283791670955125739
304 + mx2 * (0.37612638903183752464
305 + mx2 * (0.11283791670955125739
306 + mx2 * (0.026866170645131251760
307 + mx2 * 0.0052239776254421878422))));
308 #endif
309 }
310
311 // compute the error function erf(z)
FADDEEVA(erf)312 cmplx FADDEEVA(erf)(cmplx z, double relerr)
313 {
314 double x = creal(z), y = cimag(z);
315
316 if (y == 0)
317 return C(FADDEEVA_RE(erf)(x),
318 y); // preserve sign of 0
319 if (x == 0) // handle separately for speed & handling of y = Inf or NaN
320 return C(x, // preserve sign of 0
321 /* handle y -> Inf limit manually, since
322 exp(y^2) -> Inf but Im[w(y)] -> 0, so
323 IEEE will give us a NaN when it should be Inf */
324 y*y > 720 ? (y > 0 ? Inf : -Inf)
325 : exp(y*y) * FADDEEVA(w_im)(y));
326
327 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
328 double mIm_z2 = -2*x*y; // Im(-z^2)
329 if (mRe_z2 < -750) // underflow
330 return (x >= 0 ? 1.0 : -1.0);
331
332 /* Handle positive and negative x via different formulas,
333 using the mirror symmetries of w, to avoid overflow/underflow
334 problems from multiplying exponentially large and small quantities. */
335 if (x >= 0) {
336 if (x < 8e-2) {
337 if (fabs(y) < 1e-2)
338 goto taylor;
339 else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
340 goto taylor_erfi;
341 }
342 /* don't use complex exp function, since that will produce spurious NaN
343 values when multiplying w in an overflow situation. */
344 return 1.0 - exp(mRe_z2) *
345 (C(cos(mIm_z2), sin(mIm_z2))
346 * FADDEEVA(w)(C(-y,x), relerr));
347 }
348 else { // x < 0
349 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
350 if (fabs(y) < 1e-2)
351 goto taylor;
352 else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
353 goto taylor_erfi;
354 }
355 else if (isnan(x))
356 return C(NaN, y == 0 ? 0 : NaN);
357 /* don't use complex exp function, since that will produce spurious NaN
358 values when multiplying w in an overflow situation. */
359 return exp(mRe_z2) *
360 (C(cos(mIm_z2), sin(mIm_z2))
361 * FADDEEVA(w)(C(y,-x), relerr)) - 1.0;
362 }
363
364 // Use Taylor series for small |z|, to avoid cancellation inaccuracy
365 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
366 taylor:
367 {
368 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
369 return z * (1.1283791670955125739
370 + mz2 * (0.37612638903183752464
371 + mz2 * (0.11283791670955125739
372 + mz2 * (0.026866170645131251760
373 + mz2 * 0.0052239776254421878422))));
374 }
375
376 /* for small |x| and small |xy|,
377 use Taylor series to avoid cancellation inaccuracy:
378 erf(x+iy) = erf(iy)
379 + 2*exp(y^2)/sqrt(pi) *
380 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
381 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
382 where:
383 erf(iy) = exp(y^2) * Im[w(y)]
384 */
385 taylor_erfi:
386 {
387 double x2 = x*x, y2 = y*y;
388 double expy2 = exp(y2);
389 return C
390 (expy2 * x * (1.1283791670955125739
391 - x2 * (0.37612638903183752464
392 + 0.75225277806367504925*y2)
393 + x2*x2 * (0.11283791670955125739
394 + y2 * (0.45135166683820502956
395 + 0.15045055561273500986*y2))),
396 expy2 * (FADDEEVA(w_im)(y)
397 - x2*y * (1.1283791670955125739
398 - x2 * (0.56418958354775628695
399 + 0.37612638903183752464*y2))));
400 }
401 }
402
403 /////////////////////////////////////////////////////////////////////////
404
405 // return sinc(x) = sin(x)/x, given both x and sin(x)
406 // [since we only use this in cases where sin(x) has already been computed]
sinc(double x,double sinx)407 static inline double sinc(double x, double sinx) {
408 return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x;
409 }
410
411 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2
sinh_taylor(double x)412 static inline double sinh_taylor(double x) {
413 return x * (1 + (x*x) * (0.1666666666666666666667
414 + 0.00833333333333333333333 * (x*x)));
415 }
416
sqr(double x)417 static inline double sqr(double x) { return x*x; }
418
419 // precomputed table of expa2n2[n-1] = exp(-a2*n*n)
420 // for double-precision a2 = 0.26865... in FADDEEVA(w), below.
421 static const double expa2n2[] = {
422 7.64405281671221563e-01,
423 3.41424527166548425e-01,
424 8.91072646929412548e-02,
425 1.35887299055460086e-02,
426 1.21085455253437481e-03,
427 6.30452613933449404e-05,
428 1.91805156577114683e-06,
429 3.40969447714832381e-08,
430 3.54175089099469393e-10,
431 2.14965079583260682e-12,
432 7.62368911833724354e-15,
433 1.57982797110681093e-17,
434 1.91294189103582677e-20,
435 1.35344656764205340e-23,
436 5.59535712428588720e-27,
437 1.35164257972401769e-30,
438 1.90784582843501167e-34,
439 1.57351920291442930e-38,
440 7.58312432328032845e-43,
441 2.13536275438697082e-47,
442 3.51352063787195769e-52,
443 3.37800830266396920e-57,
444 1.89769439468301000e-62,
445 6.22929926072668851e-68,
446 1.19481172006938722e-73,
447 1.33908181133005953e-79,
448 8.76924303483223939e-86,
449 3.35555576166254986e-92,
450 7.50264110688173024e-99,
451 9.80192200745410268e-106,
452 7.48265412822268959e-113,
453 3.33770122566809425e-120,
454 8.69934598159861140e-128,
455 1.32486951484088852e-135,
456 1.17898144201315253e-143,
457 6.13039120236180012e-152,
458 1.86258785950822098e-160,
459 3.30668408201432783e-169,
460 3.43017280887946235e-178,
461 2.07915397775808219e-187,
462 7.36384545323984966e-197,
463 1.52394760394085741e-206,
464 1.84281935046532100e-216,
465 1.30209553802992923e-226,
466 5.37588903521080531e-237,
467 1.29689584599763145e-247,
468 1.82813078022866562e-258,
469 1.50576355348684241e-269,
470 7.24692320799294194e-281,
471 2.03797051314726829e-292,
472 3.34880215927873807e-304,
473 0.0 // underflow (also prevents reads past array end, below)
474 };
475
476 /////////////////////////////////////////////////////////////////////////
477
FADDEEVA(w)478 cmplx FADDEEVA(w)(cmplx z, double relerr)
479 {
480 if (creal(z) == 0.0)
481 return C(FADDEEVA_RE(erfcx)(cimag(z)),
482 creal(z)); // give correct sign of 0 in cimag(w)
483 else if (cimag(z) == 0)
484 return C(exp(-sqr(creal(z))),
485 FADDEEVA(w_im)(creal(z)));
486
487 double a, a2, c;
488 if (relerr <= DBL_EPSILON) {
489 relerr = DBL_EPSILON;
490 a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))
491 c = 0.329973702884629072537; // (2/pi) * a;
492 a2 = 0.268657157075235951582; // a^2
493 }
494 else {
495 const double pi = 3.14159265358979323846264338327950288419716939937510582;
496 if (relerr > 0.1) relerr = 0.1; // not sensible to compute < 1 digit
497 a = pi / sqrt(-log(relerr*0.5));
498 c = (2/pi)*a;
499 a2 = a*a;
500 }
501 const double x = fabs(creal(z));
502 const double y = cimag(z), ya = fabs(y);
503
504 cmplx ret = 0.; // return value
505
506 double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;
507
508 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z|
509
510 #if USE_CONTINUED_FRACTION
511 if (ya > 7 || (x > 6 // continued fraction is faster
512 /* As pointed out by M. Zaghloul, the continued
513 fraction seems to give a large relative error in
514 Re w(z) for |x| ~ 6 and small |y|, so use
515 algorithm 816 in this region: */
516 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
517
518 /* Poppe & Wijers suggest using a number of terms
519 nu = 3 + 1442 / (26*rho + 77)
520 where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
521 (They only use this expansion for rho >= 1, but rho a little less
522 than 1 seems okay too.)
523 Instead, I did my own fit to a slightly different function
524 that avoids the hypotenuse calculation, using NLopt to minimize
525 the sum of the squares of the errors in nu with the constraint
526 that the estimated nu be >= minimum nu to attain machine precision.
527 I also separate the regions where nu == 2 and nu == 1. */
528 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
529 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
530 if (x + ya > 4000) { // nu <= 2
531 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z
532 // scale to avoid overflow
533 if (x > ya) {
534 double yax = ya / xs;
535 double denom = ispi / (xs + yax*ya);
536 ret = C(denom*yax, denom);
537 }
538 else if (isinf(ya))
539 return ((isnan(x) || y < 0)
540 ? C(NaN,NaN) : C(0,0));
541 else {
542 double xya = xs / ya;
543 double denom = ispi / (xya*xs + ya);
544 ret = C(denom, denom*xya);
545 }
546 }
547 else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
548 double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
549 double denom = ispi / (dr*dr + di*di);
550 ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
551 }
552 }
553 else { // compute nu(z) estimate and do general continued fraction
554 const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit
555 double nu = floor(c0 + c1 / (c2*x + c3*ya + c4));
556 double wr = xs, wi = ya;
557 for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) {
558 // w <- z - nu/w:
559 double denom = nu / (wr*wr + wi*wi);
560 wr = xs - wr * denom;
561 wi = ya + wi * denom;
562 }
563 { // w(z) = i/sqrt(pi) / w:
564 double denom = ispi / (wr*wr + wi*wi);
565 ret = C(denom*wi, denom*wr);
566 }
567 }
568 if (y < 0) {
569 // use w(z) = 2.0*exp(-z*z) - w(-z),
570 // but be careful of overflow in exp(-z*z)
571 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
572 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
573 }
574 else
575 return ret;
576 }
577 #else // !USE_CONTINUED_FRACTION
578 if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision
579 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
580 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
581 // scale to avoid overflow
582 if (x > ya) {
583 double yax = ya / xs;
584 double denom = ispi / (xs + yax*ya);
585 ret = C(denom*yax, denom);
586 }
587 else {
588 double xya = xs / ya;
589 double denom = ispi / (xya*xs + ya);
590 ret = C(denom, denom*xya);
591 }
592 if (y < 0) {
593 // use w(z) = 2.0*exp(-z*z) - w(-z),
594 // but be careful of overflow in exp(-z*z)
595 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
596 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
597 }
598 else
599 return ret;
600 }
601 #endif // !USE_CONTINUED_FRACTION
602
603 /* Note: The test that seems to be suggested in the paper is x <
604 sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2)
605 underflows to zero and sum1,sum2,sum4 are zero. However, long
606 before this occurs, the sum1,sum2,sum4 contributions are
607 negligible in double precision; I find that this happens for x >
608 about 6, for all y. On the other hand, I find that the case
609 where we compute all of the sums is faster (at least with the
610 precomputed expa2n2 table) until about x=10. Furthermore, if we
611 try to compute all of the sums for x > 20, I find that we
612 sometimes run into numerical problems because underflow/overflow
613 problems start to appear in the various coefficients of the sums,
614 below. Therefore, we use x < 10 here. */
615 else if (x < 10) {
616 double prod2ax = 1, prodm2ax = 1;
617 double expx2;
618
619 if (isnan(y))
620 return C(y,y);
621
622 /* Somewhat ugly copy-and-paste duplication here, but I see significant
623 speedups from using the special-case code with the precomputed
624 exponential, and the x < 5e-4 special case is needed for accuracy. */
625
626 if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table
627 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
628 const double x2 = x*x;
629 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
630 // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
631 const double ax2 = 1.036642960860171859744*x; // 2*a*x
632 const double exp2ax =
633 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
634 const double expm2ax =
635 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
636 for (int n = 1; 1; ++n) {
637 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
638 prod2ax *= exp2ax;
639 prodm2ax *= expm2ax;
640 sum1 += coef;
641 sum2 += coef * prodm2ax;
642 sum3 += coef * prod2ax;
643
644 // really = sum5 - sum4
645 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
646
647 // test convergence via sum3
648 if (coef * prod2ax < relerr * sum3) break;
649 }
650 }
651 else { // x > 5e-4, compute sum4 and sum5 separately
652 expx2 = exp(-x*x);
653 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
654 for (int n = 1; 1; ++n) {
655 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
656 prod2ax *= exp2ax;
657 prodm2ax *= expm2ax;
658 sum1 += coef;
659 sum2 += coef * prodm2ax;
660 sum4 += (coef * prodm2ax) * (a*n);
661 sum3 += coef * prod2ax;
662 sum5 += (coef * prod2ax) * (a*n);
663 // test convergence via sum5, since this sum has the slowest decay
664 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
665 }
666 }
667 }
668 else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly
669 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
670 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
671 const double x2 = x*x;
672 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
673 for (int n = 1; 1; ++n) {
674 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
675 prod2ax *= exp2ax;
676 prodm2ax *= expm2ax;
677 sum1 += coef;
678 sum2 += coef * prodm2ax;
679 sum3 += coef * prod2ax;
680
681 // really = sum5 - sum4
682 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
683
684 // test convergence via sum3
685 if (coef * prod2ax < relerr * sum3) break;
686 }
687 }
688 else { // x > 5e-4, compute sum4 and sum5 separately
689 expx2 = exp(-x*x);
690 for (int n = 1; 1; ++n) {
691 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
692 prod2ax *= exp2ax;
693 prodm2ax *= expm2ax;
694 sum1 += coef;
695 sum2 += coef * prodm2ax;
696 sum4 += (coef * prodm2ax) * (a*n);
697 sum3 += coef * prod2ax;
698 sum5 += (coef * prod2ax) * (a*n);
699 // test convergence via sum5, since this sum has the slowest decay
700 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
701 }
702 }
703 }
704 const double expx2erfcxy = // avoid spurious overflow for large negative y
705 y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision
706 ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x);
707 if (y > 5) { // imaginary terms cancel
708 const double sinxy = sin(x*y);
709 ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
710 + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
711 }
712 else {
713 double xs = creal(z);
714 const double sinxy = sin(xs*y);
715 const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y);
716 const double coef1 = expx2erfcxy - c*y*sum1;
717 const double coef2 = c*xs*expx2;
718 ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
719 coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
720 }
721 }
722 else { // x large: only sum3 & sum5 contribute (see above note)
723 if (isnan(x))
724 return C(x,x);
725 if (isnan(y))
726 return C(y,y);
727
728 #if USE_CONTINUED_FRACTION
729 ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
730 #else
731 if (y < 0) {
732 /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so
733 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
734 if y*y - x*x > -36 or so. So, compute this term just in case.
735 We also need the -exp(-x*x) term to compute Re[w] accurately
736 in the case where y is very small. */
737 ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
738 }
739 else
740 ret = exp(-x*x); // not negligible in real part if y very small
741 #endif
742 // (round instead of ceil as in original paper; note that x/a > 1 here)
743 double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0
744 double dx = a*n0 - x;
745 sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y);
746 sum5 = a*n0 * sum3;
747 double exp1 = exp(4*a*dx), exp1dn = 1;
748 int dn;
749 for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms
750 double np = n0 + dn, nm = n0 - dn;
751 double tp = exp(-sqr(a*dn+dx));
752 double tm = tp * (exp1dn *= exp1); // trick to get tm from tp
753 tp /= (a2*(np*np) + y*y);
754 tm /= (a2*(nm*nm) + y*y);
755 sum3 += tp + tm;
756 sum5 += a * (np * tp + nm * tm);
757 if (a * (np * tp + nm * tm) < relerr * sum5) goto finish;
758 }
759 while (1) { // loop over n0+dn terms only (since n0-dn <= 0)
760 double np = n0 + dn++;
761 double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y);
762 sum3 += tp;
763 sum5 += a * np * tp;
764 if (a * np * tp < relerr * sum5) goto finish;
765 }
766 }
767 finish:
768 return ret + C((0.5*c)*y*(sum2+sum3),
769 (0.5*c)*copysign(sum5-sum4, creal(z)));
770 }
771
772 /////////////////////////////////////////////////////////////////////////
773
774 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by
775 Steven G. Johnson, October 2012.
776
777 This function combines a few different ideas.
778
779 First, for x > 50, it uses a continued-fraction expansion (same as
780 for the Faddeeva function, but with algebraic simplifications for z=i*x).
781
782 Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations,
783 but with two twists:
784
785 a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation,
786 inspired by a similar transformation in the octave-forge/specfun
787 erfcx by Soren Hauberg, results in much faster Chebyshev convergence
788 than other simple transformations I have examined.
789
790 b) Instead of using a single Chebyshev polynomial for the entire
791 [0,1] y interval, we break the interval up into 100 equal
792 subintervals, with a switch/lookup table, and use much lower
793 degree Chebyshev polynomials in each subinterval. This greatly
794 improves performance in my tests.
795
796 For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
797 with the usual checks for overflow etcetera.
798
799 Performance-wise, it seems to be substantially faster than either
800 the SLATEC DERFC function [or an erfcx function derived therefrom]
801 or Cody's CALERF function (from netlib.org/specfun), while
802 retaining near machine precision in accuracy. */
803
804 /* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x).
805
806 Uses a look-up table of 100 different Chebyshev polynomials
807 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
808 with the help of Maple and a little shell script. This allows
809 the Chebyshev polynomials to be of significantly lower degree (about 1/4)
810 compared to fitting the whole [0,1] interval with a single polynomial. */
erfcx_y100(double y100)811 static double erfcx_y100(double y100)
812 {
813 switch ((int) y100) {
814 case 0: {
815 double t = 2*y100 - 1;
816 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;
817 }
818 case 1: {
819 double t = 2*y100 - 3;
820 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;
821 }
822 case 2: {
823 double t = 2*y100 - 5;
824 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;
825 }
826 case 3: {
827 double t = 2*y100 - 7;
828 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;
829 }
830 case 4: {
831 double t = 2*y100 - 9;
832 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;
833 }
834 case 5: {
835 double t = 2*y100 - 11;
836 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;
837 }
838 case 6: {
839 double t = 2*y100 - 13;
840 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;
841 }
842 case 7: {
843 double t = 2*y100 - 15;
844 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;
845 }
846 case 8: {
847 double t = 2*y100 - 17;
848 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;
849 }
850 case 9: {
851 double t = 2*y100 - 19;
852 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;
853 }
854 case 10: {
855 double t = 2*y100 - 21;
856 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;
857 }
858 case 11: {
859 double t = 2*y100 - 23;
860 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;
861 }
862 case 12: {
863 double t = 2*y100 - 25;
864 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;
865 }
866 case 13: {
867 double t = 2*y100 - 27;
868 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;
869 }
870 case 14: {
871 double t = 2*y100 - 29;
872 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;
873 }
874 case 15: {
875 double t = 2*y100 - 31;
876 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;
877 }
878 case 16: {
879 double t = 2*y100 - 33;
880 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;
881 }
882 case 17: {
883 double t = 2*y100 - 35;
884 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;
885 }
886 case 18: {
887 double t = 2*y100 - 37;
888 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;
889 }
890 case 19: {
891 double t = 2*y100 - 39;
892 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;
893 }
894 case 20: {
895 double t = 2*y100 - 41;
896 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;
897 }
898 case 21: {
899 double t = 2*y100 - 43;
900 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;
901 }
902 case 22: {
903 double t = 2*y100 - 45;
904 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;
905 }
906 case 23: {
907 double t = 2*y100 - 47;
908 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;
909 }
910 case 24: {
911 double t = 2*y100 - 49;
912 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;
913 }
914 case 25: {
915 double t = 2*y100 - 51;
916 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;
917 }
918 case 26: {
919 double t = 2*y100 - 53;
920 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;
921 }
922 case 27: {
923 double t = 2*y100 - 55;
924 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;
925 }
926 case 28: {
927 double t = 2*y100 - 57;
928 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;
929 }
930 case 29: {
931 double t = 2*y100 - 59;
932 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;
933 }
934 case 30: {
935 double t = 2*y100 - 61;
936 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;
937 }
938 case 31: {
939 double t = 2*y100 - 63;
940 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;
941 }
942 case 32: {
943 double t = 2*y100 - 65;
944 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;
945 }
946 case 33: {
947 double t = 2*y100 - 67;
948 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;
949 }
950 case 34: {
951 double t = 2*y100 - 69;
952 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;
953 }
954 case 35: {
955 double t = 2*y100 - 71;
956 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;
957 }
958 case 36: {
959 double t = 2*y100 - 73;
960 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;
961 }
962 case 37: {
963 double t = 2*y100 - 75;
964 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;
965 }
966 case 38: {
967 double t = 2*y100 - 77;
968 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;
969 }
970 case 39: {
971 double t = 2*y100 - 79;
972 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;
973 }
974 case 40: {
975 double t = 2*y100 - 81;
976 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;
977 }
978 case 41: {
979 double t = 2*y100 - 83;
980 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;
981 }
982 case 42: {
983 double t = 2*y100 - 85;
984 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;
985 }
986 case 43: {
987 double t = 2*y100 - 87;
988 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;
989 }
990 case 44: {
991 double t = 2*y100 - 89;
992 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;
993 }
994 case 45: {
995 double t = 2*y100 - 91;
996 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;
997 }
998 case 46: {
999 double t = 2*y100 - 93;
1000 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;
1001 }
1002 case 47: {
1003 double t = 2*y100 - 95;
1004 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;
1005 }
1006 case 48: {
1007 double t = 2*y100 - 97;
1008 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;
1009 }
1010 case 49: {
1011 double t = 2*y100 - 99;
1012 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;
1013 }
1014 case 50: {
1015 double t = 2*y100 - 101;
1016 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;
1017 }
1018 case 51: {
1019 double t = 2*y100 - 103;
1020 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;
1021 }
1022 case 52: {
1023 double t = 2*y100 - 105;
1024 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;
1025 }
1026 case 53: {
1027 double t = 2*y100 - 107;
1028 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;
1029 }
1030 case 54: {
1031 double t = 2*y100 - 109;
1032 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;
1033 }
1034 case 55: {
1035 double t = 2*y100 - 111;
1036 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;
1037 }
1038 case 56: {
1039 double t = 2*y100 - 113;
1040 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;
1041 }
1042 case 57: {
1043 double t = 2*y100 - 115;
1044 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;
1045 }
1046 case 58: {
1047 double t = 2*y100 - 117;
1048 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;
1049 }
1050 case 59: {
1051 double t = 2*y100 - 119;
1052 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;
1053 }
1054 case 60: {
1055 double t = 2*y100 - 121;
1056 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;
1057 }
1058 case 61: {
1059 double t = 2*y100 - 123;
1060 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;
1061 }
1062 case 62: {
1063 double t = 2*y100 - 125;
1064 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;
1065 }
1066 case 63: {
1067 double t = 2*y100 - 127;
1068 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;
1069 }
1070 case 64: {
1071 double t = 2*y100 - 129;
1072 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;
1073 }
1074 case 65: {
1075 double t = 2*y100 - 131;
1076 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;
1077 }
1078 case 66: {
1079 double t = 2*y100 - 133;
1080 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;
1081 }
1082 case 67: {
1083 double t = 2*y100 - 135;
1084 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;
1085 }
1086 case 68: {
1087 double t = 2*y100 - 137;
1088 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;
1089 }
1090 case 69: {
1091 double t = 2*y100 - 139;
1092 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;
1093 }
1094 case 70: {
1095 double t = 2*y100 - 141;
1096 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;
1097 }
1098 case 71: {
1099 double t = 2*y100 - 143;
1100 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;
1101 }
1102 case 72: {
1103 double t = 2*y100 - 145;
1104 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;
1105 }
1106 case 73: {
1107 double t = 2*y100 - 147;
1108 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;
1109 }
1110 case 74: {
1111 double t = 2*y100 - 149;
1112 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;
1113 }
1114 case 75: {
1115 double t = 2*y100 - 151;
1116 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;
1117 }
1118 case 76: {
1119 double t = 2*y100 - 153;
1120 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;
1121 }
1122 case 77: {
1123 double t = 2*y100 - 155;
1124 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;
1125 }
1126 case 78: {
1127 double t = 2*y100 - 157;
1128 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;
1129 }
1130 case 79: {
1131 double t = 2*y100 - 159;
1132 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;
1133 }
1134 case 80: {
1135 double t = 2*y100 - 161;
1136 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;
1137 }
1138 case 81: {
1139 double t = 2*y100 - 163;
1140 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;
1141 }
1142 case 82: {
1143 double t = 2*y100 - 165;
1144 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;
1145 }
1146 case 83: {
1147 double t = 2*y100 - 167;
1148 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;
1149 }
1150 case 84: {
1151 double t = 2*y100 - 169;
1152 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;
1153 }
1154 case 85: {
1155 double t = 2*y100 - 171;
1156 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;
1157 }
1158 case 86: {
1159 double t = 2*y100 - 173;
1160 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;
1161 }
1162 case 87: {
1163 double t = 2*y100 - 175;
1164 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;
1165 }
1166 case 88: {
1167 double t = 2*y100 - 177;
1168 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;
1169 }
1170 case 89: {
1171 double t = 2*y100 - 179;
1172 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;
1173 }
1174 case 90: {
1175 double t = 2*y100 - 181;
1176 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;
1177 }
1178 case 91: {
1179 double t = 2*y100 - 183;
1180 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;
1181 }
1182 case 92: {
1183 double t = 2*y100 - 185;
1184 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;
1185 }
1186 case 93: {
1187 double t = 2*y100 - 187;
1188 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;
1189 }
1190 case 94: {
1191 double t = 2*y100 - 189;
1192 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;
1193 }
1194 case 95: {
1195 double t = 2*y100 - 191;
1196 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;
1197 }
1198 case 96: {
1199 double t = 2*y100 - 193;
1200 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;
1201 }
1202 case 97: {
1203 double t = 2*y100 - 195;
1204 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;
1205 }
1206 case 98: {
1207 double t = 2*y100 - 197;
1208 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;
1209 }
1210 case 99: {
1211 double t = 2*y100 - 199;
1212 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;
1213 }
1214 }
1215 // we only get here if y = 1, i.e. |x| < 4*eps, in which case
1216 // erfcx is within 1e-15 of 1..
1217 return 1.0;
1218 }
1219
FADDEEVA_RE(erfcx)1220 double FADDEEVA_RE(erfcx)(double x)
1221 {
1222 if (x >= 0) {
1223 if (x > 50) { // continued-fraction expansion is faster
1224 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1225 if (x > 5e7) // 1-term expansion, important to avoid overflow
1226 return ispi / x;
1227 /* 5-term expansion (rely on compiler for CSE), simplified from:
1228 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */
1229 return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
1230 }
1231 return erfcx_y100(400/(4+x));
1232 }
1233 else
1234 return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x)
1235 : 2*exp(x*x) - erfcx_y100(400/(4-x)));
1236 }
1237
1238 /////////////////////////////////////////////////////////////////////////
1239 /* Compute a scaled Dawson integral
1240 FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
1241 equivalent to the imaginary part w(x) for real x.
1242
1243 Uses methods similar to the erfcx calculation above: continued fractions
1244 for large |x|, a lookup table of Chebyshev polynomials for smaller |x|,
1245 and finally a Taylor expansion for |x|<0.01.
1246
1247 Steven G. Johnson, October 2012. */
1248
1249 /* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).
1250
1251 Uses a look-up table of 100 different Chebyshev polynomials
1252 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
1253 with the help of Maple and a little shell script. This allows
1254 the Chebyshev polynomials to be of significantly lower degree (about 1/30)
1255 compared to fitting the whole [0,1] interval with a single polynomial. */
w_im_y100(double y100,double x)1256 static double w_im_y100(double y100, double x) {
1257 switch ((int) y100) {
1258 case 0: {
1259 double t = 2*y100 - 1;
1260 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;
1261 }
1262 case 1: {
1263 double t = 2*y100 - 3;
1264 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;
1265 }
1266 case 2: {
1267 double t = 2*y100 - 5;
1268 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;
1269 }
1270 case 3: {
1271 double t = 2*y100 - 7;
1272 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;
1273 }
1274 case 4: {
1275 double t = 2*y100 - 9;
1276 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;
1277 }
1278 case 5: {
1279 double t = 2*y100 - 11;
1280 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;
1281 }
1282 case 6: {
1283 double t = 2*y100 - 13;
1284 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;
1285 }
1286 case 7: {
1287 double t = 2*y100 - 15;
1288 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;
1289 }
1290 case 8: {
1291 double t = 2*y100 - 17;
1292 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;
1293 }
1294 case 9: {
1295 double t = 2*y100 - 19;
1296 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;
1297 }
1298 case 10: {
1299 double t = 2*y100 - 21;
1300 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;
1301 }
1302 case 11: {
1303 double t = 2*y100 - 23;
1304 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;
1305 }
1306 case 12: {
1307 double t = 2*y100 - 25;
1308 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;
1309 }
1310 case 13: {
1311 double t = 2*y100 - 27;
1312 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;
1313 }
1314 case 14: {
1315 double t = 2*y100 - 29;
1316 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;
1317 }
1318 case 15: {
1319 double t = 2*y100 - 31;
1320 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;
1321 }
1322 case 16: {
1323 double t = 2*y100 - 33;
1324 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;
1325 }
1326 case 17: {
1327 double t = 2*y100 - 35;
1328 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;
1329 }
1330 case 18: {
1331 double t = 2*y100 - 37;
1332 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;
1333 }
1334 case 19: {
1335 double t = 2*y100 - 39;
1336 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;
1337 }
1338 case 20: {
1339 double t = 2*y100 - 41;
1340 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;
1341 }
1342 case 21: {
1343 double t = 2*y100 - 43;
1344 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;
1345 }
1346 case 22: {
1347 double t = 2*y100 - 45;
1348 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;
1349 }
1350 case 23: {
1351 double t = 2*y100 - 47;
1352 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;
1353 }
1354 case 24: {
1355 double t = 2*y100 - 49;
1356 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;
1357 }
1358 case 25: {
1359 double t = 2*y100 - 51;
1360 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;
1361 }
1362 case 26: {
1363 double t = 2*y100 - 53;
1364 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;
1365 }
1366 case 27: {
1367 double t = 2*y100 - 55;
1368 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;
1369 }
1370 case 28: {
1371 double t = 2*y100 - 57;
1372 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;
1373 }
1374 case 29: {
1375 double t = 2*y100 - 59;
1376 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;
1377 }
1378 case 30: {
1379 double t = 2*y100 - 61;
1380 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;
1381 }
1382 case 31: {
1383 double t = 2*y100 - 63;
1384 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;
1385 }
1386 case 32: {
1387 double t = 2*y100 - 65;
1388 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;
1389 }
1390 case 33: {
1391 double t = 2*y100 - 67;
1392 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;
1393 }
1394 case 34: {
1395 double t = 2*y100 - 69;
1396 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;
1397 }
1398 case 35: {
1399 double t = 2*y100 - 71;
1400 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;
1401 }
1402 case 36: {
1403 double t = 2*y100 - 73;
1404 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;
1405 }
1406 case 37: {
1407 double t = 2*y100 - 75;
1408 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;
1409 }
1410 case 38: {
1411 double t = 2*y100 - 77;
1412 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;
1413 }
1414 case 39: {
1415 double t = 2*y100 - 79;
1416 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;
1417 }
1418 case 40: {
1419 double t = 2*y100 - 81;
1420 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;
1421 }
1422 case 41: {
1423 double t = 2*y100 - 83;
1424 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;
1425 }
1426 case 42: {
1427 double t = 2*y100 - 85;
1428 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;
1429 }
1430 case 43: {
1431 double t = 2*y100 - 87;
1432 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;
1433 }
1434 case 44: {
1435 double t = 2*y100 - 89;
1436 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;
1437 }
1438 case 45: {
1439 double t = 2*y100 - 91;
1440 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;
1441 }
1442 case 46: {
1443 double t = 2*y100 - 93;
1444 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;
1445 }
1446 case 47: {
1447 double t = 2*y100 - 95;
1448 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;
1449 }
1450 case 48: {
1451 double t = 2*y100 - 97;
1452 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;
1453 }
1454 case 49: {
1455 double t = 2*y100 - 99;
1456 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;
1457 }
1458 case 50: {
1459 double t = 2*y100 - 101;
1460 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;
1461 }
1462 case 51: {
1463 double t = 2*y100 - 103;
1464 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;
1465 }
1466 case 52: {
1467 double t = 2*y100 - 105;
1468 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;
1469 }
1470 case 53: {
1471 double t = 2*y100 - 107;
1472 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;
1473 }
1474 case 54: {
1475 double t = 2*y100 - 109;
1476 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;
1477 }
1478 case 55: {
1479 double t = 2*y100 - 111;
1480 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;
1481 }
1482 case 56: {
1483 double t = 2*y100 - 113;
1484 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;
1485 }
1486 case 57: {
1487 double t = 2*y100 - 115;
1488 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;
1489 }
1490 case 58: {
1491 double t = 2*y100 - 117;
1492 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;
1493 }
1494 case 59: {
1495 double t = 2*y100 - 119;
1496 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;
1497 }
1498 case 60: {
1499 double t = 2*y100 - 121;
1500 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;
1501 }
1502 case 61: {
1503 double t = 2*y100 - 123;
1504 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;
1505 }
1506 case 62: {
1507 double t = 2*y100 - 125;
1508 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;
1509 }
1510 case 63: {
1511 double t = 2*y100 - 127;
1512 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;
1513 }
1514 case 64: {
1515 double t = 2*y100 - 129;
1516 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;
1517 }
1518 case 65: {
1519 double t = 2*y100 - 131;
1520 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;
1521 }
1522 case 66: {
1523 double t = 2*y100 - 133;
1524 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;
1525 }
1526 case 67: {
1527 double t = 2*y100 - 135;
1528 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;
1529 }
1530 case 68: {
1531 double t = 2*y100 - 137;
1532 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;
1533 }
1534 case 69: {
1535 double t = 2*y100 - 139;
1536 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;
1537 }
1538 case 70: {
1539 double t = 2*y100 - 141;
1540 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;
1541 }
1542 case 71: {
1543 double t = 2*y100 - 143;
1544 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;
1545 }
1546 case 72: {
1547 double t = 2*y100 - 145;
1548 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;
1549 }
1550 case 73: {
1551 double t = 2*y100 - 147;
1552 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;
1553 }
1554 case 74: {
1555 double t = 2*y100 - 149;
1556 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;
1557 }
1558 case 75: {
1559 double t = 2*y100 - 151;
1560 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;
1561 }
1562 case 76: {
1563 double t = 2*y100 - 153;
1564 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;
1565 }
1566 case 77: {
1567 double t = 2*y100 - 155;
1568 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;
1569 }
1570 case 78: {
1571 double t = 2*y100 - 157;
1572 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;
1573 }
1574 case 79: {
1575 double t = 2*y100 - 159;
1576 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;
1577 }
1578 case 80: {
1579 double t = 2*y100 - 161;
1580 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;
1581 }
1582 case 81: {
1583 double t = 2*y100 - 163;
1584 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;
1585 }
1586 case 82: {
1587 double t = 2*y100 - 165;
1588 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;
1589 }
1590 case 83: {
1591 double t = 2*y100 - 167;
1592 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;
1593 }
1594 case 84: {
1595 double t = 2*y100 - 169;
1596 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;
1597 }
1598 case 85: {
1599 double t = 2*y100 - 171;
1600 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;
1601 }
1602 case 86: {
1603 double t = 2*y100 - 173;
1604 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;
1605 }
1606 case 87: {
1607 double t = 2*y100 - 175;
1608 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;
1609 }
1610 case 88: {
1611 double t = 2*y100 - 177;
1612 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;
1613 }
1614 case 89: {
1615 double t = 2*y100 - 179;
1616 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;
1617 }
1618 case 90: {
1619 double t = 2*y100 - 181;
1620 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;
1621 }
1622 case 91: {
1623 double t = 2*y100 - 183;
1624 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;
1625 }
1626 case 92: {
1627 double t = 2*y100 - 185;
1628 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;
1629 }
1630 case 93: {
1631 double t = 2*y100 - 187;
1632 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;
1633 }
1634 case 94: {
1635 double t = 2*y100 - 189;
1636 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;
1637 }
1638 case 95: {
1639 double t = 2*y100 - 191;
1640 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;
1641 }
1642 case 96: {
1643 double t = 2*y100 - 193;
1644 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;
1645 }
1646 case 97: case 98:
1647 case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
1648 // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9)
1649 double x2 = x*x;
1650 return x * (1.1283791670955125739
1651 - x2 * (0.75225277806367504925
1652 - x2 * (0.30090111122547001970
1653 - x2 * (0.085971746064420005629
1654 - x2 * 0.016931216931216931217))));
1655 }
1656 }
1657 /* Since 0 <= y100 < 101, this is only reached if x is NaN,
1658 in which case we should return NaN. */
1659 return NaN;
1660 }
1661
FADDEEVA(w_im)1662 double FADDEEVA(w_im)(double x)
1663 {
1664 if (x >= 0) {
1665 if (x > 45) { // continued-fraction expansion is faster
1666 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1667 if (x > 5e7) // 1-term expansion, important to avoid overflow
1668 return ispi / x;
1669 /* 5-term expansion (rely on compiler for CSE), simplified from:
1670 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1671 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1672 }
1673 return w_im_y100(100/(1+x), x);
1674 }
1675 else { // = -FADDEEVA(w_im)(-x)
1676 if (x < -45) { // continued-fraction expansion is faster
1677 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1678 if (x < -5e7) // 1-term expansion, important to avoid overflow
1679 return ispi / x;
1680 /* 5-term expansion (rely on compiler for CSE), simplified from:
1681 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1682 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1683 }
1684 return -w_im_y100(100/(1-x), -x);
1685 }
1686 }
1687