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