1 // -*- mode:c++; tab-width:2; indent-tabs-mode:nil; -*-
2
3 /* Copyright (c) 2012 Massachusetts Institute of Technology
4 *
5 * Permission is hereby granted, free of charge, to any person obtaining
6 * a copy of this software and associated documentation files (the
7 * "Software"), to deal in the Software without restriction, including
8 * without limitation the rights to use, copy, modify, merge, publish,
9 * distribute, sublicense, and/or sell copies of the Software, and to
10 * permit persons to whom the Software is furnished to do so, subject to
11 * the following conditions:
12 *
13 * The above copyright notice and this permission notice shall be
14 * included in all copies or substantial portions of the Software.
15 *
16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
17 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
18 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
19 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
20 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
21 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
22 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
23 */
24
25 /* (Note that this file can be compiled with either C++, in which
26 case it uses C++ std::complex<double>, or C, in which case it
27 uses C99 double complex.) */
28
29 /* Available at: http://ab-initio.mit.edu/Faddeeva
30
31 Computes various error functions (erf, erfc, erfi, erfcx),
32 including the Dawson integral, in the complex plane, based
33 on algorithms for the computation of the Faddeeva function
34 w(z) = exp(-z^2) * erfc(-i*z).
35 Given w(z), the error functions are mostly straightforward
36 to compute, except for certain regions where we have to
37 switch to Taylor expansions to avoid cancellation errors
38 [e.g., near the origin for erf(z)].
39
40 To compute the Faddeeva function, we use a combination of two
41 algorithms:
42
43 For sufficiently large |z|, we use a continued-fraction expansion
44 for w(z) similar to those described in:
45
46 Walter Gautschi, "Efficient computation of the complex error
47 function," SIAM J. Numer. Anal. 7(1), pp. 187-198 (1970)
48
49 G. P. M. Poppe and C. M. J. Wijers, "More efficient computation
50 of the complex error function," ACM Trans. Math. Soft. 16(1),
51 pp. 38-46 (1990).
52
53 Unlike those papers, however, we switch to a completely different
54 algorithm for smaller |z|:
55
56 Mofreh R. Zaghloul and Ahmed N. Ali, "Algorithm 916: Computing the
57 Faddeyeva and Voigt Functions," ACM Trans. Math. Soft. 38(2), 15
58 (2011).
59
60 (I initially used this algorithm for all z, but it turned out to be
61 significantly slower than the continued-fraction expansion for
62 larger |z|. On the other hand, it is competitive for smaller |z|,
63 and is significantly more accurate than the Poppe & Wijers code
64 in some regions, e.g., in the vicinity of z=1+1i.)
65
66 Note that this is an INDEPENDENT RE-IMPLEMENTATION of these algorithms,
67 based on the description in the papers ONLY. In particular, I did
68 not refer to the authors' Fortran or Matlab implementations, respectively,
69 (which are under restrictive ACM copyright terms and therefore unusable
70 in free/open-source software).
71
72 Steven G. Johnson, Massachusetts Institute of Technology
73 http://math.mit.edu/~stevenj
74 October 2012.
75
76 -- Note that Algorithm 916 assumes that the erfc(x) function,
77 or rather the scaled function erfcx(x) = exp(x*x)*erfc(x),
78 is supplied for REAL arguments x. I originally used an
79 erfcx routine derived from DERFC in SLATEC, but I have
80 since replaced it with a much faster routine written by
81 me which uses a combination of continued-fraction expansions
82 and a lookup table of Chebyshev polynomials. For speed,
83 I implemented a similar algorithm for Im[w(x)] of real x,
84 since this comes up frequently in the other error functions.
85
86 A small test program is included the end, which checks
87 the w(z) etc. results against several known values. To compile
88 the test function, compile with -DTEST_FADDEEVA (that is,
89 #define TEST_FADDEEVA).
90
91 If HAVE_CONFIG_H is #defined (e.g., by compiling with -DHAVE_CONFIG_H),
92 then we #include "config.h", which is assumed to be a GNU autoconf-style
93 header defining HAVE_* macros to indicate the presence of features. In
94 particular, if HAVE_ISNAN and HAVE_ISINF are #defined, we use those
95 functions in math.h instead of defining our own, and if HAVE_ERF and/or
96 HAVE_ERFC are defined we use those functions from <cmath> for erf and
97 erfc of real arguments, respectively, instead of defining our own.
98
99 REVISION HISTORY:
100 4 October 2012: Initial public release (SGJ)
101 5 October 2012: Revised (SGJ) to fix spelling error,
102 start summation for large x at round(x/a) (> 1)
103 rather than ceil(x/a) as in the original
104 paper, which should slightly improve performance
105 (and, apparently, slightly improves accuracy)
106 19 October 2012: Revised (SGJ) to fix bugs for large x, large -y,
107 and 15<x<26. Performance improvements. Prototype
108 now supplies default value for relerr.
109 24 October 2012: Switch to continued-fraction expansion for
110 sufficiently large z, for performance reasons.
111 Also, avoid spurious overflow for |z| > 1e154.
112 Set relerr argument to min(relerr,0.1).
113 27 October 2012: Enhance accuracy in Re[w(z)] taken by itself,
114 by switching to Alg. 916 in a region near
115 the real-z axis where continued fractions
116 have poor relative accuracy in Re[w(z)]. Thanks
117 to M. Zaghloul for the tip.
118 29 October 2012: Replace SLATEC-derived erfcx routine with
119 completely rewritten code by me, using a very
120 different algorithm which is much faster.
121 30 October 2012: Implemented special-case code for real z
122 (where real part is exp(-x^2) and imag part is
123 Dawson integral), using algorithm similar to erfx.
124 Export ImFaddeeva_w function to make Dawson's
125 integral directly accessible.
126 3 November 2012: Provide implementations of erf, erfc, erfcx,
127 and Dawson functions in Faddeeva:: namespace,
128 in addition to Faddeeva::w. Provide header
129 file Faddeeva.hh.
130 4 November 2012: Slightly faster erf for real arguments.
131 Updated MATLAB and Octave plugins.
132 27 November 2012: Support compilation with either C++ or
133 plain C (using C99 complex numbers).
134 For real x, use standard-library erf(x)
135 and erfc(x) if available (for C99 or C++11).
136 #include "config.h" if HAVE_CONFIG_H is #defined.
137 15 December 2012: Portability fixes (copysign, Inf/NaN creation),
138 use CMPLX/__builtin_complex if available in C,
139 slight accuracy improvements to erf and dawson
140 functions near the origin. Use gnulib functions
141 if GNULIB_NAMESPACE is defined.
142 18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson)
143 */
144
145 /////////////////////////////////////////////////////////////////////////
146 /* If this file is compiled as a part of a larger project,
147 support using an autoconf-style config.h header file
148 (with various "HAVE_*" #defines to indicate features)
149 if HAVE_CONFIG_H is #defined (in GNU autotools style). */
150
151 #if defined (HAVE_CONFIG_H)
152 # include "config.h"
153 #endif
154
155 /////////////////////////////////////////////////////////////////////////
156 // macros to allow us to use either C++ or C (with C99 features)
157
158 #if defined (__cplusplus)
159
160 # include "lo-ieee.h"
161
162 # include "Faddeeva.hh"
163
164 # include <cfloat>
165 # include <cmath>
166 # include <limits>
167
168 // use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS)
169 # define Inf octave::numeric_limits<double>::Inf ()
170 # define NaN octave::numeric_limits<double>::NaN ()
171
172 typedef std::complex<double> cmplx;
173
174 // Use C-like complex syntax, since the C syntax is more restrictive
175 # define cexp(z) exp(z)
176 # define creal(z) real(z)
177 # define cimag(z) imag(z)
178 # define cpolar(r,t) polar(r,t)
179
180 # define C(a,b) cmplx(a,b)
181
182 # define FADDEEVA(name) Faddeeva::name
183 # define FADDEEVA_RE(name) Faddeeva::name
184
185 // isnan/isinf were introduced in C++11
186 # if defined (lo_ieee_isnan) && defined (lo_ieee_isinf)
187 # define isnan lo_ieee_isnan
188 # define isinf lo_ieee_isinf
189 # elif (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF))
my_isnan(double x)190 static inline bool my_isnan(double x) { return x != x; }
191 # define isnan my_isnan
my_isinf(double x)192 static inline bool my_isinf(double x) { return 1/x == 0.; }
193 # define isinf my_isinf
194 # elif (__cplusplus >= 201103L)
195 // g++ gets confused between the C and C++ isnan/isinf functions
196 # define isnan std::isnan
197 # define isinf std::isinf
198 # endif
199
200 // copysign was introduced in C++11 (and is also in POSIX and C99)
201 # if defined(_WIN32) || defined(__WIN32__)
202 # define copysign _copysign // of course MS had to be different
203 # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX)
my_copysign(double x,double y)204 static inline double my_copysign(double x, double y) { return y<0 ? -x : x; }
205 # define copysign my_copysign
206 # endif
207
208 #else // !__cplusplus, i.e., pure C (requires C99 features)
209
210 # include "Faddeeva.h"
211
212 # define _GNU_SOURCE // enable GNU libc NAN extension if possible
213
214 # include <float.h>
215 # include <math.h>
216
217 typedef double complex cmplx;
218
219 # define FADDEEVA(name) Faddeeva_ ## name
220 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
221
222 /* Constructing complex numbers like 0+i*NaN is problematic in C99
223 without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
224 I is a complex (rather than imaginary) constant. For some reason,
225 however, it works fine in (pre-4.7) gcc if I define Inf and NaN as
226 1/0 and 0/0 (and only if I compile with optimization -O1 or more),
227 but not if I use the INFINITY or NAN macros. */
228
229 /* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro
230 may not be defined unless we are using a recent (2012) version of
231 glibc and compile with -std=c11... note that icc lies about being
232 gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
233 # if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER))
234 # define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
235 # endif
236
237 # if defined (CMPLX) // C11
238 # define C(a,b) CMPLX(a,b)
239 # define Inf INFINITY // C99 infinity
240 # if defined (NAN) // GNU libc extension
241 # define NaN NAN
242 # else
243 # define NaN (0./0.) // NaN
244 # endif
245 # else
246 # define C(a,b) ((a) + I*(b))
247 # define Inf (1./0.)
248 # define NaN (0./0.)
249 # endif
250
cpolar(double r,double t)251 static inline cmplx cpolar(double r, double t)
252 {
253 if (r == 0.0 && !isnan(t))
254 return 0.0;
255 else
256 return C(r * cos(t), r * sin(t));
257 }
258
259 #endif // !__cplusplus, i.e., pure C (requires C99 features)
260
261 /////////////////////////////////////////////////////////////////////////
262 // Auxiliary routines to compute other special functions based on w(z)
263
264 // compute erfcx(z) = exp(z^2) erfz(z)
FADDEEVA(erfcx)265 cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
266 {
267 return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
268 }
269
270 // compute the error function erf(x)
FADDEEVA_RE(erf)271 double FADDEEVA_RE(erf)(double x)
272 {
273 #if !defined(__cplusplus)
274 return erf(x); // C99 supplies erf in math.h
275 #elif (__cplusplus >= 201103L) || defined(HAVE_ERF)
276 return ::erf(x); // C++11 supplies std::erf in cmath
277 #else
278 double mx2 = -x*x;
279 if (mx2 < -750) // underflow
280 return (x >= 0 ? 1.0 : -1.0);
281
282 if (x >= 0) {
283 if (x < 8e-2) goto taylor;
284 return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
285 }
286 else { // x < 0
287 if (x > -8e-2) goto taylor;
288 return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0;
289 }
290
291 // Use Taylor series for small |x|, to avoid cancellation inaccuracy
292 // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...)
293 taylor:
294 return x * (1.1283791670955125739
295 + mx2 * (0.37612638903183752464
296 + mx2 * (0.11283791670955125739
297 + mx2 * (0.026866170645131251760
298 + mx2 * 0.0052239776254421878422))));
299 #endif
300 }
301
302 // compute the error function erf(z)
FADDEEVA(erf)303 cmplx FADDEEVA(erf)(cmplx z, double relerr)
304 {
305 double x = creal(z), y = cimag(z);
306
307 if (y == 0)
308 return C(FADDEEVA_RE(erf)(x),
309 y); // preserve sign of 0
310 if (x == 0) // handle separately for speed & handling of y = Inf or NaN
311 return C(x, // preserve sign of 0
312 /* handle y -> Inf limit manually, since
313 exp(y^2) -> Inf but Im[w(y)] -> 0, so
314 IEEE will give us a NaN when it should be Inf */
315 y*y > 720 ? (y > 0 ? Inf : -Inf)
316 : exp(y*y) * FADDEEVA(w_im)(y));
317
318 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
319 double mIm_z2 = -2*x*y; // Im(-z^2)
320 if (mRe_z2 < -750) // underflow
321 return (x >= 0 ? 1.0 : -1.0);
322
323 /* Handle positive and negative x via different formulas,
324 using the mirror symmetries of w, to avoid overflow/underflow
325 problems from multiplying exponentially large and small quantities. */
326 if (x >= 0) {
327 if (x < 8e-2) {
328 if (fabs(y) < 1e-2)
329 goto taylor;
330 else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
331 goto taylor_erfi;
332 }
333 /* don't use complex exp function, since that will produce spurious NaN
334 values when multiplying w in an overflow situation. */
335 return 1.0 - exp(mRe_z2) *
336 (C(cos(mIm_z2), sin(mIm_z2))
337 * FADDEEVA(w)(C(-y,x), relerr));
338 }
339 else { // x < 0
340 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
341 if (fabs(y) < 1e-2)
342 goto taylor;
343 else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
344 goto taylor_erfi;
345 }
346 else if (isnan(x))
347 return C(NaN, y == 0 ? 0 : NaN);
348 /* don't use complex exp function, since that will produce spurious NaN
349 values when multiplying w in an overflow situation. */
350 return exp(mRe_z2) *
351 (C(cos(mIm_z2), sin(mIm_z2))
352 * FADDEEVA(w)(C(y,-x), relerr)) - 1.0;
353 }
354
355 // Use Taylor series for small |z|, to avoid cancellation inaccuracy
356 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
357 taylor:
358 {
359 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
360 return z * (1.1283791670955125739
361 + mz2 * (0.37612638903183752464
362 + mz2 * (0.11283791670955125739
363 + mz2 * (0.026866170645131251760
364 + mz2 * 0.0052239776254421878422))));
365 }
366
367 /* for small |x| and small |xy|,
368 use Taylor series to avoid cancellation inaccuracy:
369 erf(x+iy) = erf(iy)
370 + 2*exp(y^2)/sqrt(pi) *
371 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
372 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
373 where:
374 erf(iy) = exp(y^2) * Im[w(y)]
375 */
376 taylor_erfi:
377 {
378 double x2 = x*x, y2 = y*y;
379 double expy2 = exp(y2);
380 return C
381 (expy2 * x * (1.1283791670955125739
382 - x2 * (0.37612638903183752464
383 + 0.75225277806367504925*y2)
384 + x2*x2 * (0.11283791670955125739
385 + y2 * (0.45135166683820502956
386 + 0.15045055561273500986*y2))),
387 expy2 * (FADDEEVA(w_im)(y)
388 - x2*y * (1.1283791670955125739
389 - x2 * (0.56418958354775628695
390 + 0.37612638903183752464*y2))));
391 }
392 }
393
394 // erfi(z) = -i erf(iz)
FADDEEVA(erfi)395 cmplx FADDEEVA(erfi)(cmplx z, double relerr)
396 {
397 cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr);
398 return C(cimag(e), -creal(e));
399 }
400
401 // erfi(x) = -i erf(ix)
FADDEEVA_RE(erfi)402 double FADDEEVA_RE(erfi)(double x)
403 {
404 return x*x > 720 ? (x > 0 ? Inf : -Inf)
405 : exp(x*x) * FADDEEVA(w_im)(x);
406 }
407
408 // erfc(x) = 1 - erf(x)
FADDEEVA_RE(erfc)409 double FADDEEVA_RE(erfc)(double x)
410 {
411 #if !defined(__cplusplus)
412 return erfc(x); // C99 supplies erfc in math.h
413 #elif (__cplusplus >= 201103L) || defined(HAVE_ERFC)
414 return ::erfc(x); // C++11 supplies std::erfc in cmath
415 #else
416 if (x*x > 750) // underflow
417 return (x >= 0 ? 0.0 : 2.0);
418 return (x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x)
419 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x));
420 #endif
421 }
422
423 // erfc(z) = 1 - erf(z)
FADDEEVA(erfc)424 cmplx FADDEEVA(erfc)(cmplx z, double relerr)
425 {
426 double x = creal(z), y = cimag(z);
427
428 if (x == 0.)
429 return C(1,
430 /* handle y -> Inf limit manually, since
431 exp(y^2) -> Inf but Im[w(y)] -> 0, so
432 IEEE will give us a NaN when it should be Inf */
433 y*y > 720 ? (y > 0 ? -Inf : Inf)
434 : -exp(y*y) * FADDEEVA(w_im)(y));
435 if (y == 0.) {
436 if (x*x > 750) // underflow
437 return C(x >= 0 ? 0.0 : 2.0,
438 -y); // preserve sign of 0
439 return C(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x)
440 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x),
441 -y); // preserve sign of zero
442 }
443
444 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
445 double mIm_z2 = -2*x*y; // Im(-z^2)
446 if (mRe_z2 < -750) // underflow
447 return (x >= 0 ? 0.0 : 2.0);
448
449 if (x >= 0)
450 return cexp(C(mRe_z2, mIm_z2))
451 * FADDEEVA(w)(C(-y,x), relerr);
452 else
453 return 2.0 - cexp(C(mRe_z2, mIm_z2))
454 * FADDEEVA(w)(C(y,-x), relerr);
455 }
456
457 // compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x)
FADDEEVA_RE(Dawson)458 double FADDEEVA_RE(Dawson)(double x)
459 {
460 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
461 return spi2 * FADDEEVA(w_im)(x);
462 }
463
464 // compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z)
FADDEEVA(Dawson)465 cmplx FADDEEVA(Dawson)(cmplx z, double relerr)
466 {
467 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
468 double x = creal(z), y = cimag(z);
469
470 // handle axes separately for speed & proper handling of x or y = Inf or NaN
471 if (y == 0)
472 return C(spi2 * FADDEEVA(w_im)(x),
473 -y); // preserve sign of 0
474 if (x == 0) {
475 double y2 = y*y;
476 if (y2 < 2.5e-5) { // Taylor expansion
477 return C(x, // preserve sign of 0
478 y * (1.
479 + y2 * (0.6666666666666666666666666666666666666667
480 + y2 * 0.26666666666666666666666666666666666667)));
481 }
482 return C(x, // preserve sign of 0
483 spi2 * (y >= 0
484 ? exp(y2) - FADDEEVA_RE(erfcx)(y)
485 : FADDEEVA_RE(erfcx)(-y) - exp(y2)));
486 }
487
488 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
489 double mIm_z2 = -2*x*y; // Im(-z^2)
490 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
491
492 /* Handle positive and negative x via different formulas,
493 using the mirror symmetries of w, to avoid overflow/underflow
494 problems from multiplying exponentially large and small quantities. */
495 if (y >= 0) {
496 if (y < 5e-3) {
497 if (fabs(x) < 5e-3)
498 goto taylor;
499 else if (fabs(mIm_z2) < 5e-3)
500 goto taylor_realaxis;
501 }
502 cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr);
503 return spi2 * C(-cimag(res), creal(res));
504 }
505 else { // y < 0
506 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
507 if (fabs(x) < 5e-3)
508 goto taylor;
509 else if (fabs(mIm_z2) < 5e-3)
510 goto taylor_realaxis;
511 }
512 else if (isnan(y))
513 return C(x == 0 ? 0 : NaN, NaN);
514 cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2);
515 return spi2 * C(-cimag(res), creal(res));
516 }
517
518 // Use Taylor series for small |z|, to avoid cancellation inaccuracy
519 // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
520 taylor:
521 return z * (1.
522 + mz2 * (0.6666666666666666666666666666666666666667
523 + mz2 * 0.2666666666666666666666666666666666666667));
524
525 /* for small |y| and small |xy|,
526 use Taylor series to avoid cancellation inaccuracy:
527 dawson(x + iy)
528 = D + y^2 (D + x - 2Dx^2)
529 + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3)
530 + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3)
531 + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ...
532 where D = dawson(x)
533
534 However, for large |x|, 2Dx -> 1 which gives cancellation problems in
535 this series (many of the leading terms cancel). So, for large |x|,
536 we need to substitute a continued-fraction expansion for D.
537
538 dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...))))))
539
540 The 6 terms shown here seems to be the minimum needed to be
541 accurate as soon as the simpler Taylor expansion above starts
542 breaking down. Using this 6-term expansion, factoring out the
543 denominator, and simplifying with Maple, we obtain:
544
545 Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x
546 = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4
547 Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y
548 = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4
549
550 Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction
551 expansion for the real part, and a 2-term expansion for the imaginary
552 part. (This avoids overflow problems for huge |x|.) This yields:
553
554 Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x)
555 Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1)
556
557 */
558 taylor_realaxis:
559 {
560 double x2 = x*x;
561 if (x2 > 1600) { // |x| > 40
562 double y2 = y*y;
563 if (x2 > 25e14) {// |x| > 5e7
564 double xy2 = (x*y)*(x*y);
565 return C((0.5 + y2 * (0.5 + 0.25*y2
566 - 0.16666666666666666667*xy2)) / x,
567 y * (-1 + y2 * (-0.66666666666666666667
568 + 0.13333333333333333333*xy2
569 - 0.26666666666666666667*y2))
570 / (2*x2 - 1));
571 }
572 return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
573 C(x * (33 + x2 * (-28 + 4*x2)
574 + y2 * (18 - 4*x2 + 4*y2)),
575 y * (-15 + x2 * (24 - 4*x2)
576 + y2 * (4*x2 - 10 - 4*y2)));
577 }
578 else {
579 double D = spi2 * FADDEEVA(w_im)(x);
580 double y2 = y*y;
581 return C
582 (D + y2 * (D + x - 2*D*x2)
583 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
584 + x * (0.83333333333333333333
585 - 0.33333333333333333333 * x2)),
586 y * (1 - 2*D*x
587 + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
588 + y2*y2 * (0.26666666666666666667 -
589 x2 * (0.6 - 0.13333333333333333333 * x2)
590 - D*x * (1 - x2 * (1.3333333333333333333
591 - 0.26666666666666666667 * x2)))));
592 }
593 }
594 }
595
596 /////////////////////////////////////////////////////////////////////////
597
598 // return sinc(x) = sin(x)/x, given both x and sin(x)
599 // [since we only use this in cases where sin(x) has already been computed]
sinc(double x,double sinx)600 static inline double sinc(double x, double sinx) {
601 return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x;
602 }
603
604 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2
sinh_taylor(double x)605 static inline double sinh_taylor(double x) {
606 return x * (1 + (x*x) * (0.1666666666666666666667
607 + 0.00833333333333333333333 * (x*x)));
608 }
609
sqr(double x)610 static inline double sqr(double x) { return x*x; }
611
612 // precomputed table of expa2n2[n-1] = exp(-a2*n*n)
613 // for double-precision a2 = 0.26865... in FADDEEVA(w), below.
614 static const double expa2n2[] = {
615 7.64405281671221563e-01,
616 3.41424527166548425e-01,
617 8.91072646929412548e-02,
618 1.35887299055460086e-02,
619 1.21085455253437481e-03,
620 6.30452613933449404e-05,
621 1.91805156577114683e-06,
622 3.40969447714832381e-08,
623 3.54175089099469393e-10,
624 2.14965079583260682e-12,
625 7.62368911833724354e-15,
626 1.57982797110681093e-17,
627 1.91294189103582677e-20,
628 1.35344656764205340e-23,
629 5.59535712428588720e-27,
630 1.35164257972401769e-30,
631 1.90784582843501167e-34,
632 1.57351920291442930e-38,
633 7.58312432328032845e-43,
634 2.13536275438697082e-47,
635 3.51352063787195769e-52,
636 3.37800830266396920e-57,
637 1.89769439468301000e-62,
638 6.22929926072668851e-68,
639 1.19481172006938722e-73,
640 1.33908181133005953e-79,
641 8.76924303483223939e-86,
642 3.35555576166254986e-92,
643 7.50264110688173024e-99,
644 9.80192200745410268e-106,
645 7.48265412822268959e-113,
646 3.33770122566809425e-120,
647 8.69934598159861140e-128,
648 1.32486951484088852e-135,
649 1.17898144201315253e-143,
650 6.13039120236180012e-152,
651 1.86258785950822098e-160,
652 3.30668408201432783e-169,
653 3.43017280887946235e-178,
654 2.07915397775808219e-187,
655 7.36384545323984966e-197,
656 1.52394760394085741e-206,
657 1.84281935046532100e-216,
658 1.30209553802992923e-226,
659 5.37588903521080531e-237,
660 1.29689584599763145e-247,
661 1.82813078022866562e-258,
662 1.50576355348684241e-269,
663 7.24692320799294194e-281,
664 2.03797051314726829e-292,
665 3.34880215927873807e-304,
666 0.0 // underflow (also prevents reads past array end, below)
667 };
668
669 /////////////////////////////////////////////////////////////////////////
670
FADDEEVA(w)671 cmplx FADDEEVA(w)(cmplx z, double relerr)
672 {
673 if (creal(z) == 0.0)
674 return C(FADDEEVA_RE(erfcx)(cimag(z)),
675 creal(z)); // give correct sign of 0 in cimag(w)
676 else if (cimag(z) == 0)
677 return C(exp(-sqr(creal(z))),
678 FADDEEVA(w_im)(creal(z)));
679
680 double a, a2, c;
681 if (relerr <= DBL_EPSILON) {
682 relerr = DBL_EPSILON;
683 a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))
684 c = 0.329973702884629072537; // (2/pi) * a;
685 a2 = 0.268657157075235951582; // a^2
686 }
687 else {
688 const double pi = 3.14159265358979323846264338327950288419716939937510582;
689 if (relerr > 0.1) relerr = 0.1; // not sensible to compute < 1 digit
690 a = pi / sqrt(-log(relerr*0.5));
691 c = (2/pi)*a;
692 a2 = a*a;
693 }
694 const double x = fabs(creal(z));
695 const double y = cimag(z), ya = fabs(y);
696
697 cmplx ret = 0.; // return value
698
699 double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0;
700
701 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z|
702
703 #if USE_CONTINUED_FRACTION
704 if (ya > 7 || (x > 6 // continued fraction is faster
705 /* As pointed out by M. Zaghloul, the continued
706 fraction seems to give a large relative error in
707 Re w(z) for |x| ~ 6 and small |y|, so use
708 algorithm 816 in this region: */
709 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
710
711 /* Poppe & Wijers suggest using a number of terms
712 nu = 3 + 1442 / (26*rho + 77)
713 where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
714 (They only use this expansion for rho >= 1, but rho a little less
715 than 1 seems okay too.)
716 Instead, I did my own fit to a slightly different function
717 that avoids the hypotenuse calculation, using NLopt to minimize
718 the sum of the squares of the errors in nu with the constraint
719 that the estimated nu be >= minimum nu to attain machine precision.
720 I also separate the regions where nu == 2 and nu == 1. */
721 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
722 double xs = (y < 0 ? -creal(z) : creal(z)); // compute for -z if y < 0
723 if (x + ya > 4000) { // nu <= 2
724 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z
725 // scale to avoid overflow
726 if (x > ya) {
727 double yax = ya / xs;
728 double denom = ispi / (xs + yax*ya);
729 ret = C(denom*yax, denom);
730 }
731 else if (isinf(ya))
732 return ((isnan(x) || y < 0)
733 ? C(NaN,NaN) : C(0,0));
734 else {
735 double xya = xs / ya;
736 double denom = ispi / (xya*xs + ya);
737 ret = C(denom, denom*xya);
738 }
739 }
740 else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
741 double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
742 double denom = ispi / (dr*dr + di*di);
743 ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
744 }
745 }
746 else { // compute nu(z) estimate and do general continued fraction
747 const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit
748 double nu = floor(c0 + c1 / (c2*x + c3*ya + c4));
749 double wr = xs, wi = ya;
750 for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) {
751 // w <- z - nu/w:
752 double denom = nu / (wr*wr + wi*wi);
753 wr = xs - wr * denom;
754 wi = ya + wi * denom;
755 }
756 { // w(z) = i/sqrt(pi) / w:
757 double denom = ispi / (wr*wr + wi*wi);
758 ret = C(denom*wi, denom*wr);
759 }
760 }
761 if (y < 0) {
762 // use w(z) = 2.0*exp(-z*z) - w(-z),
763 // but be careful of overflow in exp(-z*z)
764 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
765 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
766 }
767 else
768 return ret;
769 }
770 #else // !USE_CONTINUED_FRACTION
771 if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision
772 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
773 double xs = (y < 0 ? -creal(z) : creal(z)); // compute for -z if y < 0
774 // scale to avoid overflow
775 if (x > ya) {
776 double yax = ya / xs;
777 double denom = ispi / (xs + yax*ya);
778 ret = C(denom*yax, denom);
779 }
780 else {
781 double xya = xs / ya;
782 double denom = ispi / (xya*xs + ya);
783 ret = C(denom, denom*xya);
784 }
785 if (y < 0) {
786 // use w(z) = 2.0*exp(-z*z) - w(-z),
787 // but be careful of overflow in exp(-z*z)
788 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
789 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
790 }
791 else
792 return ret;
793 }
794 #endif // !USE_CONTINUED_FRACTION
795
796 /* Note: The test that seems to be suggested in the paper is x <
797 sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2)
798 underflows to zero and sum1,sum2,sum4 are zero. However, long
799 before this occurs, the sum1,sum2,sum4 contributions are
800 negligible in double precision; I find that this happens for x >
801 about 6, for all y. On the other hand, I find that the case
802 where we compute all of the sums is faster (at least with the
803 precomputed expa2n2 table) until about x=10. Furthermore, if we
804 try to compute all of the sums for x > 20, I find that we
805 sometimes run into numerical problems because underflow/overflow
806 problems start to appear in the various coefficients of the sums,
807 below. Therefore, we use x < 10 here. */
808 else if (x < 10) {
809 double prod2ax = 1, prodm2ax = 1;
810 double expx2;
811
812 if (isnan(y))
813 return C(y,y);
814
815 /* Somewhat ugly copy-and-paste duplication here, but I see significant
816 speedups from using the special-case code with the precomputed
817 exponential, and the x < 5e-4 special case is needed for accuracy. */
818
819 if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table
820 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
821 const double x2 = x*x;
822 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
823 // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
824 const double ax2 = 1.036642960860171859744*x; // 2*a*x
825 const double exp2ax =
826 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
827 const double expm2ax =
828 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
829 for (int n = 1; 1; ++n) {
830 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
831 prod2ax *= exp2ax;
832 prodm2ax *= expm2ax;
833 sum1 += coef;
834 sum2 += coef * prodm2ax;
835 sum3 += coef * prod2ax;
836
837 // really = sum5 - sum4
838 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
839
840 // test convergence via sum3
841 if (coef * prod2ax < relerr * sum3) break;
842 }
843 }
844 else { // x > 5e-4, compute sum4 and sum5 separately
845 expx2 = exp(-x*x);
846 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
847 for (int n = 1; 1; ++n) {
848 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
849 prod2ax *= exp2ax;
850 prodm2ax *= expm2ax;
851 sum1 += coef;
852 sum2 += coef * prodm2ax;
853 sum4 += (coef * prodm2ax) * (a*n);
854 sum3 += coef * prod2ax;
855 sum5 += (coef * prod2ax) * (a*n);
856 // test convergence via sum5, since this sum has the slowest decay
857 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
858 }
859 }
860 }
861 else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly
862 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
863 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
864 const double x2 = x*x;
865 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
866 for (int n = 1; 1; ++n) {
867 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
868 prod2ax *= exp2ax;
869 prodm2ax *= expm2ax;
870 sum1 += coef;
871 sum2 += coef * prodm2ax;
872 sum3 += coef * prod2ax;
873
874 // really = sum5 - sum4
875 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
876
877 // test convergence via sum3
878 if (coef * prod2ax < relerr * sum3) break;
879 }
880 }
881 else { // x > 5e-4, compute sum4 and sum5 separately
882 expx2 = exp(-x*x);
883 for (int n = 1; 1; ++n) {
884 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
885 prod2ax *= exp2ax;
886 prodm2ax *= expm2ax;
887 sum1 += coef;
888 sum2 += coef * prodm2ax;
889 sum4 += (coef * prodm2ax) * (a*n);
890 sum3 += coef * prod2ax;
891 sum5 += (coef * prod2ax) * (a*n);
892 // test convergence via sum5, since this sum has the slowest decay
893 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
894 }
895 }
896 }
897 const double expx2erfcxy = // avoid spurious overflow for large negative y
898 y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision
899 ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x);
900 if (y > 5) { // imaginary terms cancel
901 const double sinxy = sin(x*y);
902 ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
903 + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
904 }
905 else {
906 double xs = creal(z);
907 const double sinxy = sin(xs*y);
908 const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y);
909 const double coef1 = expx2erfcxy - c*y*sum1;
910 const double coef2 = c*xs*expx2;
911 ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
912 coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
913 }
914 }
915 else { // x large: only sum3 & sum5 contribute (see above note)
916 if (isnan(x))
917 return C(x,x);
918 if (isnan(y))
919 return C(y,y);
920
921 #if USE_CONTINUED_FRACTION
922 ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
923 #else
924 if (y < 0) {
925 /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so
926 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
927 if y*y - x*x > -36 or so. So, compute this term just in case.
928 We also need the -exp(-x*x) term to compute Re[w] accurately
929 in the case where y is very small. */
930 ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
931 }
932 else
933 ret = exp(-x*x); // not negligible in real part if y very small
934 #endif
935 // (round instead of ceil as in original paper; note that x/a > 1 here)
936 double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0
937 double dx = a*n0 - x;
938 sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y);
939 sum5 = a*n0 * sum3;
940 double exp1 = exp(4*a*dx), exp1dn = 1;
941 int dn;
942 for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms
943 double np = n0 + dn, nm = n0 - dn;
944 double tp = exp(-sqr(a*dn+dx));
945 double tm = tp * (exp1dn *= exp1); // trick to get tm from tp
946 tp /= (a2*(np*np) + y*y);
947 tm /= (a2*(nm*nm) + y*y);
948 sum3 += tp + tm;
949 sum5 += a * (np * tp + nm * tm);
950 if (a * (np * tp + nm * tm) < relerr * sum5) goto finish;
951 }
952 while (1) { // loop over n0+dn terms only (since n0-dn <= 0)
953 double np = n0 + dn++;
954 double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y);
955 sum3 += tp;
956 sum5 += a * np * tp;
957 if (a * np * tp < relerr * sum5) goto finish;
958 }
959 }
960 finish:
961 return ret + C((0.5*c)*y*(sum2+sum3),
962 (0.5*c)*copysign(sum5-sum4, creal(z)));
963 }
964
965 /////////////////////////////////////////////////////////////////////////
966
967 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by
968 Steven G. Johnson, October 2012.
969
970 This function combines a few different ideas.
971
972 First, for x > 50, it uses a continued-fraction expansion (same as
973 for the Faddeeva function, but with algebraic simplifications for z=i*x).
974
975 Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations,
976 but with two twists:
977
978 a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation,
979 inspired by a similar transformation in the octave-forge/specfun
980 erfcx by Soren Hauberg, results in much faster Chebyshev convergence
981 than other simple transformations I have examined.
982
983 b) Instead of using a single Chebyshev polynomial for the entire
984 [0,1] y interval, we break the interval up into 100 equal
985 subintervals, with a switch/lookup table, and use much lower
986 degree Chebyshev polynomials in each subinterval. This greatly
987 improves performance in my tests.
988
989 For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
990 with the usual checks for overflow etcetera.
991
992 Performance-wise, it seems to be substantially faster than either
993 the SLATEC DERFC function [or an erfcx function derived therefrom]
994 or Cody's CALERF function (from netlib.org/specfun), while
995 retaining near machine precision in accuracy. */
996
997 /* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x).
998
999 Uses a look-up table of 100 different Chebyshev polynomials
1000 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
1001 with the help of Maple and a little shell script. This allows
1002 the Chebyshev polynomials to be of significantly lower degree (about 1/4)
1003 compared to fitting the whole [0,1] interval with a single polynomial. */
erfcx_y100(double y100)1004 static double erfcx_y100(double y100)
1005 {
1006 switch (static_cast<int> (y100)) {
1007 case 0: {
1008 double t = 2*y100 - 1;
1009 return 0.70878032454106438663e-3 + (0.71234091047026302958e-3 + (0.35779077297597742384e-5 + (0.17403143962587937815e-7 + (0.81710660047307788845e-10 + (0.36885022360434957634e-12 + 0.15917038551111111111e-14 * t) * t) * t) * t) * t) * t;
1010 }
1011 case 1: {
1012 double t = 2*y100 - 3;
1013 return 0.21479143208285144230e-2 + (0.72686402367379996033e-3 + (0.36843175430938995552e-5 + (0.18071841272149201685e-7 + (0.85496449296040325555e-10 + (0.38852037518534291510e-12 + 0.16868473576888888889e-14 * t) * t) * t) * t) * t) * t;
1014 }
1015 case 2: {
1016 double t = 2*y100 - 5;
1017 return 0.36165255935630175090e-2 + (0.74182092323555510862e-3 + (0.37948319957528242260e-5 + (0.18771627021793087350e-7 + (0.89484715122415089123e-10 + (0.40935858517772440862e-12 + 0.17872061464888888889e-14 * t) * t) * t) * t) * t) * t;
1018 }
1019 case 3: {
1020 double t = 2*y100 - 7;
1021 return 0.51154983860031979264e-2 + (0.75722840734791660540e-3 + (0.39096425726735703941e-5 + (0.19504168704300468210e-7 + (0.93687503063178993915e-10 + (0.43143925959079664747e-12 + 0.18939926435555555556e-14 * t) * t) * t) * t) * t) * t;
1022 }
1023 case 4: {
1024 double t = 2*y100 - 9;
1025 return 0.66457513172673049824e-2 + (0.77310406054447454920e-3 + (0.40289510589399439385e-5 + (0.20271233238288381092e-7 + (0.98117631321709100264e-10 + (0.45484207406017752971e-12 + 0.20076352213333333333e-14 * t) * t) * t) * t) * t) * t;
1026 }
1027 case 5: {
1028 double t = 2*y100 - 11;
1029 return 0.82082389970241207883e-2 + (0.78946629611881710721e-3 + (0.41529701552622656574e-5 + (0.21074693344544655714e-7 + (0.10278874108587317989e-9 + (0.47965201390613339638e-12 + 0.21285907413333333333e-14 * t) * t) * t) * t) * t) * t;
1030 }
1031 case 6: {
1032 double t = 2*y100 - 13;
1033 return 0.98039537275352193165e-2 + (0.80633440108342840956e-3 + (0.42819241329736982942e-5 + (0.21916534346907168612e-7 + (0.10771535136565470914e-9 + (0.50595972623692822410e-12 + 0.22573462684444444444e-14 * t) * t) * t) * t) * t) * t;
1034 }
1035 case 7: {
1036 double t = 2*y100 - 15;
1037 return 0.11433927298290302370e-1 + (0.82372858383196561209e-3 + (0.44160495311765438816e-5 + (0.22798861426211986056e-7 + (0.11291291745879239736e-9 + (0.53386189365816880454e-12 + 0.23944209546666666667e-14 * t) * t) * t) * t) * t) * t;
1038 }
1039 case 8: {
1040 double t = 2*y100 - 17;
1041 return 0.13099232878814653979e-1 + (0.84167002467906968214e-3 + (0.45555958988457506002e-5 + (0.23723907357214175198e-7 + (0.11839789326602695603e-9 + (0.56346163067550237877e-12 + 0.25403679644444444444e-14 * t) * t) * t) * t) * t) * t;
1042 }
1043 case 9: {
1044 double t = 2*y100 - 19;
1045 return 0.14800987015587535621e-1 + (0.86018092946345943214e-3 + (0.47008265848816866105e-5 + (0.24694040760197315333e-7 + (0.12418779768752299093e-9 + (0.59486890370320261949e-12 + 0.26957764568888888889e-14 * t) * t) * t) * t) * t) * t;
1046 }
1047 case 10: {
1048 double t = 2*y100 - 21;
1049 return 0.16540351739394069380e-1 + (0.87928458641241463952e-3 + (0.48520195793001753903e-5 + (0.25711774900881709176e-7 + (0.13030128534230822419e-9 + (0.62820097586874779402e-12 + 0.28612737351111111111e-14 * t) * t) * t) * t) * t) * t;
1050 }
1051 case 11: {
1052 double t = 2*y100 - 23;
1053 return 0.18318536789842392647e-1 + (0.89900542647891721692e-3 + (0.50094684089553365810e-5 + (0.26779777074218070482e-7 + (0.13675822186304615566e-9 + (0.66358287745352705725e-12 + 0.30375273884444444444e-14 * t) * t) * t) * t) * t) * t;
1054 }
1055 case 12: {
1056 double t = 2*y100 - 25;
1057 return 0.20136801964214276775e-1 + (0.91936908737673676012e-3 + (0.51734830914104276820e-5 + (0.27900878609710432673e-7 + (0.14357976402809042257e-9 + (0.70114790311043728387e-12 + 0.32252476000000000000e-14 * t) * t) * t) * t) * t) * t;
1058 }
1059 case 13: {
1060 double t = 2*y100 - 27;
1061 return 0.21996459598282740954e-1 + (0.94040248155366777784e-3 + (0.53443911508041164739e-5 + (0.29078085538049374673e-7 + (0.15078844500329731137e-9 + (0.74103813647499204269e-12 + 0.34251892320000000000e-14 * t) * t) * t) * t) * t) * t;
1062 }
1063 case 14: {
1064 double t = 2*y100 - 29;
1065 return 0.23898877187226319502e-1 + (0.96213386835900177540e-3 + (0.55225386998049012752e-5 + (0.30314589961047687059e-7 + (0.15840826497296335264e-9 + (0.78340500472414454395e-12 + 0.36381553564444444445e-14 * t) * t) * t) * t) * t) * t;
1066 }
1067 case 15: {
1068 double t = 2*y100 - 31;
1069 return 0.25845480155298518485e-1 + (0.98459293067820123389e-3 + (0.57082915920051843672e-5 + (0.31613782169164830118e-7 + (0.16646478745529630813e-9 + (0.82840985928785407942e-12 + 0.38649975768888888890e-14 * t) * t) * t) * t) * t) * t;
1070 }
1071 case 16: {
1072 double t = 2*y100 - 33;
1073 return 0.27837754783474696598e-1 + (0.10078108563256892757e-2 + (0.59020366493792212221e-5 + (0.32979263553246520417e-7 + (0.17498524159268458073e-9 + (0.87622459124842525110e-12 + 0.41066206488888888890e-14 * t) * t) * t) * t) * t) * t;
1074 }
1075 case 17: {
1076 double t = 2*y100 - 35;
1077 return 0.29877251304899307550e-1 + (0.10318204245057349310e-2 + (0.61041829697162055093e-5 + (0.34414860359542720579e-7 + (0.18399863072934089607e-9 + (0.92703227366365046533e-12 + 0.43639844053333333334e-14 * t) * t) * t) * t) * t) * t;
1078 }
1079 case 18: {
1080 double t = 2*y100 - 37;
1081 return 0.31965587178596443475e-1 + (0.10566560976716574401e-2 + (0.63151633192414586770e-5 + (0.35924638339521924242e-7 + (0.19353584758781174038e-9 + (0.98102783859889264382e-12 + 0.46381060817777777779e-14 * t) * t) * t) * t) * t) * t;
1082 }
1083 case 19: {
1084 double t = 2*y100 - 39;
1085 return 0.34104450552588334840e-1 + (0.10823541191350532574e-2 + (0.65354356159553934436e-5 + (0.37512918348533521149e-7 + (0.20362979635817883229e-9 + (0.10384187833037282363e-11 + 0.49300625262222222221e-14 * t) * t) * t) * t) * t) * t;
1086 }
1087 case 20: {
1088 double t = 2*y100 - 41;
1089 return 0.36295603928292425716e-1 + (0.11089526167995268200e-2 + (0.67654845095518363577e-5 + (0.39184292949913591646e-7 + (0.21431552202133775150e-9 + (0.10994259106646731797e-11 + 0.52409949102222222221e-14 * t) * t) * t) * t) * t) * t;
1090 }
1091 case 21: {
1092 double t = 2*y100 - 43;
1093 return 0.38540888038840509795e-1 + (0.11364917134175420009e-2 + (0.70058230641246312003e-5 + (0.40943644083718586939e-7 + (0.22563034723692881631e-9 + (0.11642841011361992885e-11 + 0.55721092871111111110e-14 * t) * t) * t) * t) * t) * t;
1094 }
1095 case 22: {
1096 double t = 2*y100 - 45;
1097 return 0.40842225954785960651e-1 + (0.11650136437945673891e-2 + (0.72569945502343006619e-5 + (0.42796161861855042273e-7 + (0.23761401711005024162e-9 + (0.12332431172381557035e-11 + 0.59246802364444444445e-14 * t) * t) * t) * t) * t) * t;
1098 }
1099 case 23: {
1100 double t = 2*y100 - 47;
1101 return 0.43201627431540222422e-1 + (0.11945628793917272199e-2 + (0.75195743532849206263e-5 + (0.44747364553960993492e-7 + (0.25030885216472953674e-9 + (0.13065684400300476484e-11 + 0.63000532853333333334e-14 * t) * t) * t) * t) * t) * t;
1102 }
1103 case 24: {
1104 double t = 2*y100 - 49;
1105 return 0.45621193513810471438e-1 + (0.12251862608067529503e-2 + (0.77941720055551920319e-5 + (0.46803119830954460212e-7 + (0.26375990983978426273e-9 + (0.13845421370977119765e-11 + 0.66996477404444444445e-14 * t) * t) * t) * t) * t) * t;
1106 }
1107 case 25: {
1108 double t = 2*y100 - 51;
1109 return 0.48103121413299865517e-1 + (0.12569331386432195113e-2 + (0.80814333496367673980e-5 + (0.48969667335682018324e-7 + (0.27801515481905748484e-9 + (0.14674637611609884208e-11 + 0.71249589351111111110e-14 * t) * t) * t) * t) * t) * t;
1110 }
1111 case 26: {
1112 double t = 2*y100 - 53;
1113 return 0.50649709676983338501e-1 + (0.12898555233099055810e-2 + (0.83820428414568799654e-5 + (0.51253642652551838659e-7 + (0.29312563849675507232e-9 + (0.15556512782814827846e-11 + 0.75775607822222222221e-14 * t) * t) * t) * t) * t) * t;
1114 }
1115 case 27: {
1116 double t = 2*y100 - 55;
1117 return 0.53263363664388864181e-1 + (0.13240082443256975769e-2 + (0.86967260015007658418e-5 + (0.53662102750396795566e-7 + (0.30914568786634796807e-9 + (0.16494420240828493176e-11 + 0.80591079644444444445e-14 * t) * t) * t) * t) * t) * t;
1118 }
1119 case 28: {
1120 double t = 2*y100 - 57;
1121 return 0.55946601353500013794e-1 + (0.13594491197408190706e-2 + (0.90262520233016380987e-5 + (0.56202552975056695376e-7 + (0.32613310410503135996e-9 + (0.17491936862246367398e-11 + 0.85713381688888888890e-14 * t) * t) * t) * t) * t) * t;
1122 }
1123 case 29: {
1124 double t = 2*y100 - 59;
1125 return 0.58702059496154081813e-1 + (0.13962391363223647892e-2 + (0.93714365487312784270e-5 + (0.58882975670265286526e-7 + (0.34414937110591753387e-9 + (0.18552853109751857859e-11 + 0.91160736711111111110e-14 * t) * t) * t) * t) * t) * t;
1126 }
1127 case 30: {
1128 double t = 2*y100 - 61;
1129 return 0.61532500145144778048e-1 + (0.14344426411912015247e-2 + (0.97331446201016809696e-5 + (0.61711860507347175097e-7 + (0.36325987418295300221e-9 + (0.19681183310134518232e-11 + 0.96952238400000000000e-14 * t) * t) * t) * t) * t) * t;
1130 }
1131 case 31: {
1132 double t = 2*y100 - 63;
1133 return 0.64440817576653297993e-1 + (0.14741275456383131151e-2 + (0.10112293819576437838e-4 + (0.64698236605933246196e-7 + (0.38353412915303665586e-9 + (0.20881176114385120186e-11 + 0.10310784480000000000e-13 * t) * t) * t) * t) * t) * t;
1134 }
1135 case 32: {
1136 double t = 2*y100 - 65;
1137 return 0.67430045633130393282e-1 + (0.15153655418916540370e-2 + (0.10509857606888328667e-4 + (0.67851706529363332855e-7 + (0.40504602194811140006e-9 + (0.22157325110542534469e-11 + 0.10964842115555555556e-13 * t) * t) * t) * t) * t) * t;
1138 }
1139 case 33: {
1140 double t = 2*y100 - 67;
1141 return 0.70503365513338850709e-1 + (0.15582323336495709827e-2 + (0.10926868866865231089e-4 + (0.71182482239613507542e-7 + (0.42787405890153386710e-9 + (0.23514379522274416437e-11 + 0.11659571751111111111e-13 * t) * t) * t) * t) * t) * t;
1142 }
1143 case 34: {
1144 double t = 2*y100 - 69;
1145 return 0.73664114037944596353e-1 + (0.16028078812438820413e-2 + (0.11364423678778207991e-4 + (0.74701423097423182009e-7 + (0.45210162777476488324e-9 + (0.24957355004088569134e-11 + 0.12397238257777777778e-13 * t) * t) * t) * t) * t) * t;
1146 }
1147 case 35: {
1148 double t = 2*y100 - 71;
1149 return 0.76915792420819562379e-1 + (0.16491766623447889354e-2 + (0.11823685320041302169e-4 + (0.78420075993781544386e-7 + (0.47781726956916478925e-9 + (0.26491544403815724749e-11 + 0.13180196462222222222e-13 * t) * t) * t) * t) * t) * t;
1150 }
1151 case 36: {
1152 double t = 2*y100 - 73;
1153 return 0.80262075578094612819e-1 + (0.16974279491709504117e-2 + (0.12305888517309891674e-4 + (0.82350717698979042290e-7 + (0.50511496109857113929e-9 + (0.28122528497626897696e-11 + 0.14010889635555555556e-13 * t) * t) * t) * t) * t) * t;
1154 }
1155 case 37: {
1156 double t = 2*y100 - 75;
1157 return 0.83706822008980357446e-1 + (0.17476561032212656962e-2 + (0.12812343958540763368e-4 + (0.86506399515036435592e-7 + (0.53409440823869467453e-9 + (0.29856186620887555043e-11 + 0.14891851591111111111e-13 * t) * t) * t) * t) * t) * t;
1158 }
1159 case 38: {
1160 double t = 2*y100 - 77;
1161 return 0.87254084284461718231e-1 + (0.17999608886001962327e-2 + (0.13344443080089492218e-4 + (0.90900994316429008631e-7 + (0.56486134972616465316e-9 + (0.31698707080033956934e-11 + 0.15825697795555555556e-13 * t) * t) * t) * t) * t) * t;
1162 }
1163 case 39: {
1164 double t = 2*y100 - 79;
1165 return 0.90908120182172748487e-1 + (0.18544478050657699758e-2 + (0.13903663143426120077e-4 + (0.95549246062549906177e-7 + (0.59752787125242054315e-9 + (0.33656597366099099413e-11 + 0.16815130613333333333e-13 * t) * t) * t) * t) * t) * t;
1166 }
1167 case 40: {
1168 double t = 2*y100 - 81;
1169 return 0.94673404508075481121e-1 + (0.19112284419887303347e-2 + (0.14491572616545004930e-4 + (0.10046682186333613697e-6 + (0.63221272959791000515e-9 + (0.35736693975589130818e-11 + 0.17862931591111111111e-13 * t) * t) * t) * t) * t) * t;
1170 }
1171 case 41: {
1172 double t = 2*y100 - 83;
1173 return 0.98554641648004456555e-1 + (0.19704208544725622126e-2 + (0.15109836875625443935e-4 + (0.10567036667675984067e-6 + (0.66904168640019354565e-9 + (0.37946171850824333014e-11 + 0.18971959040000000000e-13 * t) * t) * t) * t) * t) * t;
1174 }
1175 case 42: {
1176 double t = 2*y100 - 85;
1177 return 0.10255677889470089531e0 + (0.20321499629472857418e-2 + (0.15760224242962179564e-4 + (0.11117756071353507391e-6 + (0.70814785110097658502e-9 + (0.40292553276632563925e-11 + 0.20145143075555555556e-13 * t) * t) * t) * t) * t) * t;
1178 }
1179 case 43: {
1180 double t = 2*y100 - 87;
1181 return 0.10668502059865093318e0 + (0.20965479776148731610e-2 + (0.16444612377624983565e-4 + (0.11700717962026152749e-6 + (0.74967203250938418991e-9 + (0.42783716186085922176e-11 + 0.21385479360000000000e-13 * t) * t) * t) * t) * t) * t;
1182 }
1183 case 44: {
1184 double t = 2*y100 - 89;
1185 return 0.11094484319386444474e0 + (0.21637548491908170841e-2 + (0.17164995035719657111e-4 + (0.12317915750735938089e-6 + (0.79376309831499633734e-9 + (0.45427901763106353914e-11 + 0.22696025653333333333e-13 * t) * t) * t) * t) * t) * t;
1186 }
1187 case 45: {
1188 double t = 2*y100 - 91;
1189 return 0.11534201115268804714e0 + (0.22339187474546420375e-2 + (0.17923489217504226813e-4 + (0.12971465288245997681e-6 + (0.84057834180389073587e-9 + (0.48233721206418027227e-11 + 0.24079890062222222222e-13 * t) * t) * t) * t) * t) * t;
1190 }
1191 case 46: {
1192 double t = 2*y100 - 93;
1193 return 0.11988259392684094740e0 + (0.23071965691918689601e-2 + (0.18722342718958935446e-4 + (0.13663611754337957520e-6 + (0.89028385488493287005e-9 + (0.51210161569225846701e-11 + 0.25540227111111111111e-13 * t) * t) * t) * t) * t) * t;
1194 }
1195 case 47: {
1196 double t = 2*y100 - 95;
1197 return 0.12457298393509812907e0 + (0.23837544771809575380e-2 + (0.19563942105711612475e-4 + (0.14396736847739470782e-6 + (0.94305490646459247016e-9 + (0.54366590583134218096e-11 + 0.27080225920000000000e-13 * t) * t) * t) * t) * t) * t;
1198 }
1199 case 48: {
1200 double t = 2*y100 - 97;
1201 return 0.12941991566142438816e0 + (0.24637684719508859484e-2 + (0.20450821127475879816e-4 + (0.15173366280523906622e-6 + (0.99907632506389027739e-9 + (0.57712760311351625221e-11 + 0.28703099555555555556e-13 * t) * t) * t) * t) * t) * t;
1202 }
1203 case 49: {
1204 double t = 2*y100 - 99;
1205 return 0.13443048593088696613e0 + (0.25474249981080823877e-2 + (0.21385669591362915223e-4 + (0.15996177579900443030e-6 + (0.10585428844575134013e-8 + (0.61258809536787882989e-11 + 0.30412080142222222222e-13 * t) * t) * t) * t) * t) * t;
1206 }
1207 case 50: {
1208 double t = 2*y100 - 101;
1209 return 0.13961217543434561353e0 + (0.26349215871051761416e-2 + (0.22371342712572567744e-4 + (0.16868008199296822247e-6 + (0.11216596910444996246e-8 + (0.65015264753090890662e-11 + 0.32210394506666666666e-13 * t) * t) * t) * t) * t) * t;
1210 }
1211 case 51: {
1212 double t = 2*y100 - 103;
1213 return 0.14497287157673800690e0 + (0.27264675383982439814e-2 + (0.23410870961050950197e-4 + (0.17791863939526376477e-6 + (0.11886425714330958106e-8 + (0.68993039665054288034e-11 + 0.34101266222222222221e-13 * t) * t) * t) * t) * t) * t;
1214 }
1215 case 52: {
1216 double t = 2*y100 - 105;
1217 return 0.15052089272774618151e0 + (0.28222846410136238008e-2 + (0.24507470422713397006e-4 + (0.18770927679626136909e-6 + (0.12597184587583370712e-8 + (0.73203433049229821618e-11 + 0.36087889048888888890e-13 * t) * t) * t) * t) * t) * t;
1218 }
1219 case 53: {
1220 double t = 2*y100 - 107;
1221 return 0.15626501395774612325e0 + (0.29226079376196624949e-2 + (0.25664553693768450545e-4 + (0.19808568415654461964e-6 + (0.13351257759815557897e-8 + (0.77658124891046760667e-11 + 0.38173420035555555555e-13 * t) * t) * t) * t) * t) * t;
1222 }
1223 case 54: {
1224 double t = 2*y100 - 109;
1225 return 0.16221449434620737567e0 + (0.30276865332726475672e-2 + (0.26885741326534564336e-4 + (0.20908350604346384143e-6 + (0.14151148144240728728e-8 + (0.82369170665974313027e-11 + 0.40360957457777777779e-13 * t) * t) * t) * t) * t) * t;
1226 }
1227 case 55: {
1228 double t = 2*y100 - 111;
1229 return 0.16837910595412130659e0 + (0.31377844510793082301e-2 + (0.28174873844911175026e-4 + (0.22074043807045782387e-6 + (0.14999481055996090039e-8 + (0.87348993661930809254e-11 + 0.42653528977777777779e-13 * t) * t) * t) * t) * t) * t;
1230 }
1231 case 56: {
1232 double t = 2*y100 - 113;
1233 return 0.17476916455659369953e0 + (0.32531815370903068316e-2 + (0.29536024347344364074e-4 + (0.23309632627767074202e-6 + (0.15899007843582444846e-8 + (0.92610375235427359475e-11 + 0.45054073102222222221e-13 * t) * t) * t) * t) * t) * t;
1234 }
1235 case 57: {
1236 double t = 2*y100 - 115;
1237 return 0.18139556223643701364e0 + (0.33741744168096996041e-2 + (0.30973511714709500836e-4 + (0.24619326937592290996e-6 + (0.16852609412267750744e-8 + (0.98166442942854895573e-11 + 0.47565418097777777779e-13 * t) * t) * t) * t) * t) * t;
1238 }
1239 case 58: {
1240 double t = 2*y100 - 117;
1241 return 0.18826980194443664549e0 + (0.35010775057740317997e-2 + (0.32491914440014267480e-4 + (0.26007572375886319028e-6 + (0.17863299617388376116e-8 + (0.10403065638343878679e-10 + 0.50190265831111111110e-13 * t) * t) * t) * t) * t) * t;
1242 }
1243 case 59: {
1244 double t = 2*y100 - 119;
1245 return 0.19540403413693967350e0 + (0.36342240767211326315e-2 + (0.34096085096200907289e-4 + (0.27479061117017637474e-6 + (0.18934228504790032826e-8 + (0.11021679075323598664e-10 + 0.52931171733333333334e-13 * t) * t) * t) * t) * t) * t;
1246 }
1247 case 60: {
1248 double t = 2*y100 - 121;
1249 return 0.20281109560651886959e0 + (0.37739673859323597060e-2 + (0.35791165457592409054e-4 + (0.29038742889416172404e-6 + (0.20068685374849001770e-8 + (0.11673891799578381999e-10 + 0.55790523093333333334e-13 * t) * t) * t) * t) * t) * t;
1250 }
1251 case 61: {
1252 double t = 2*y100 - 123;
1253 return 0.21050455062669334978e0 + (0.39206818613925652425e-2 + (0.37582602289680101704e-4 + (0.30691836231886877385e-6 + (0.21270101645763677824e-8 + (0.12361138551062899455e-10 + 0.58770520160000000000e-13 * t) * t) * t) * t) * t) * t;
1254 }
1255 case 62: {
1256 double t = 2*y100 - 125;
1257 return 0.21849873453703332479e0 + (0.40747643554689586041e-2 + (0.39476163820986711501e-4 + (0.32443839970139918836e-6 + (0.22542053491518680200e-8 + (0.13084879235290858490e-10 + 0.61873153262222222221e-13 * t) * t) * t) * t) * t) * t;
1258 }
1259 case 63: {
1260 double t = 2*y100 - 127;
1261 return 0.22680879990043229327e0 + (0.42366354648628516935e-2 + (0.41477956909656896779e-4 + (0.34300544894502810002e-6 + (0.23888264229264067658e-8 + (0.13846596292818514601e-10 + 0.65100183751111111110e-13 * t) * t) * t) * t) * t) * t;
1262 }
1263 case 64: {
1264 double t = 2*y100 - 129;
1265 return 0.23545076536988703937e0 + (0.44067409206365170888e-2 + (0.43594444916224700881e-4 + (0.36268045617760415178e-6 + (0.25312606430853202748e-8 + (0.14647791812837903061e-10 + 0.68453122631111111110e-13 * t) * t) * t) * t) * t) * t;
1266 }
1267 case 65: {
1268 double t = 2*y100 - 131;
1269 return 0.24444156740777432838e0 + (0.45855530511605787178e-2 + (0.45832466292683085475e-4 + (0.38352752590033030472e-6 + (0.26819103733055603460e-8 + (0.15489984390884756993e-10 + 0.71933206364444444445e-13 * t) * t) * t) * t) * t) * t;
1270 }
1271 case 66: {
1272 double t = 2*y100 - 133;
1273 return 0.25379911500634264643e0 + (0.47735723208650032167e-2 + (0.48199253896534185372e-4 + (0.40561404245564732314e-6 + (0.28411932320871165585e-8 + (0.16374705736458320149e-10 + 0.75541379822222222221e-13 * t) * t) * t) * t) * t) * t;
1274 }
1275 case 67: {
1276 double t = 2*y100 - 135;
1277 return 0.26354234756393613032e0 + (0.49713289477083781266e-2 + (0.50702455036930367504e-4 + (0.42901079254268185722e-6 + (0.30095422058900481753e-8 + (0.17303497025347342498e-10 + 0.79278273368888888890e-13 * t) * t) * t) * t) * t) * t;
1278 }
1279 case 68: {
1280 double t = 2*y100 - 137;
1281 return 0.27369129607732343398e0 + (0.51793846023052643767e-2 + (0.53350152258326602629e-4 + (0.45379208848865015485e-6 + (0.31874057245814381257e-8 + (0.18277905010245111046e-10 + 0.83144182364444444445e-13 * t) * t) * t) * t) * t) * t;
1282 }
1283 case 69: {
1284 double t = 2*y100 - 139;
1285 return 0.28426714781640316172e0 + (0.53983341916695141966e-2 + (0.56150884865255810638e-4 + (0.48003589196494734238e-6 + (0.33752476967570796349e-8 + (0.19299477888083469086e-10 + 0.87139049137777777779e-13 * t) * t) * t) * t) * t) * t;
1286 }
1287 case 70: {
1288 double t = 2*y100 - 141;
1289 return 0.29529231465348519920e0 + (0.56288077305420795663e-2 + (0.59113671189913307427e-4 + (0.50782393781744840482e-6 + (0.35735475025851713168e-8 + (0.20369760937017070382e-10 + 0.91262442613333333334e-13 * t) * t) * t) * t) * t) * t;
1290 }
1291 case 71: {
1292 double t = 2*y100 - 143;
1293 return 0.30679050522528838613e0 + (0.58714723032745403331e-2 + (0.62248031602197686791e-4 + (0.53724185766200945789e-6 + (0.37827999418960232678e-8 + (0.21490291930444538307e-10 + 0.95513539182222222221e-13 * t) * t) * t) * t) * t) * t;
1294 }
1295 case 72: {
1296 double t = 2*y100 - 145;
1297 return 0.31878680111173319425e0 + (0.61270341192339103514e-2 + (0.65564012259707640976e-4 + (0.56837930287837738996e-6 + (0.40035151353392378882e-8 + (0.22662596341239294792e-10 + 0.99891109760000000000e-13 * t) * t) * t) * t) * t) * t;
1298 }
1299 case 73: {
1300 double t = 2*y100 - 147;
1301 return 0.33130773722152622027e0 + (0.63962406646798080903e-2 + (0.69072209592942396666e-4 + (0.60133006661885941812e-6 + (0.42362183765883466691e-8 + (0.23888182347073698382e-10 + 0.10439349811555555556e-12 * t) * t) * t) * t) * t) * t;
1302 }
1303 case 74: {
1304 double t = 2*y100 - 149;
1305 return 0.34438138658041336523e0 + (0.66798829540414007258e-2 + (0.72783795518603561144e-4 + (0.63619220443228800680e-6 + (0.44814499336514453364e-8 + (0.25168535651285475274e-10 + 0.10901861383111111111e-12 * t) * t) * t) * t) * t) * t;
1306 }
1307 case 75: {
1308 double t = 2*y100 - 151;
1309 return 0.35803744972380175583e0 + (0.69787978834882685031e-2 + (0.76710543371454822497e-4 + (0.67306815308917386747e-6 + (0.47397647975845228205e-8 + (0.26505114141143050509e-10 + 0.11376390933333333333e-12 * t) * t) * t) * t) * t) * t;
1310 }
1311 case 76: {
1312 double t = 2*y100 - 153;
1313 return 0.37230734890119724188e0 + (0.72938706896461381003e-2 + (0.80864854542670714092e-4 + (0.71206484718062688779e-6 + (0.50117323769745883805e-8 + (0.27899342394100074165e-10 + 0.11862637614222222222e-12 * t) * t) * t) * t) * t) * t;
1314 }
1315 case 77: {
1316 double t = 2*y100 - 155;
1317 return 0.38722432730555448223e0 + (0.76260375162549802745e-2 + (0.85259785810004603848e-4 + (0.75329383305171327677e-6 + (0.52979361368388119355e-8 + (0.29352606054164086709e-10 + 0.12360253370666666667e-12 * t) * t) * t) * t) * t) * t;
1318 }
1319 case 78: {
1320 double t = 2*y100 - 157;
1321 return 0.40282355354616940667e0 + (0.79762880915029728079e-2 + (0.89909077342438246452e-4 + (0.79687137961956194579e-6 + (0.55989731807360403195e-8 + (0.30866246101464869050e-10 + 0.12868841946666666667e-12 * t) * t) * t) * t) * t) * t;
1322 }
1323 case 79: {
1324 double t = 2*y100 - 159;
1325 return 0.41914223158913787649e0 + (0.83456685186950463538e-2 + (0.94827181359250161335e-4 + (0.84291858561783141014e-6 + (0.59154537751083485684e-8 + (0.32441553034347469291e-10 + 0.13387957943111111111e-12 * t) * t) * t) * t) * t) * t;
1326 }
1327 case 80: {
1328 double t = 2*y100 - 161;
1329 return 0.43621971639463786896e0 + (0.87352841828289495773e-2 + (0.10002929142066799966e-3 + (0.89156148280219880024e-6 + (0.62480008150788597147e-8 + (0.34079760983458878910e-10 + 0.13917107176888888889e-12 * t) * t) * t) * t) * t) * t;
1330 }
1331 case 81: {
1332 double t = 2*y100 - 163;
1333 return 0.45409763548534330981e0 + (0.91463027755548240654e-2 + (0.10553137232446167258e-3 + (0.94293113464638623798e-6 + (0.65972492312219959885e-8 + (0.35782041795476563662e-10 + 0.14455745872000000000e-12 * t) * t) * t) * t) * t) * t;
1334 }
1335 case 82: {
1336 double t = 2*y100 - 165;
1337 return 0.47282001668512331468e0 + (0.95799574408860463394e-2 + (0.11135019058000067469e-3 + (0.99716373005509038080e-6 + (0.69638453369956970347e-8 + (0.37549499088161345850e-10 + 0.15003280712888888889e-12 * t) * t) * t) * t) * t) * t;
1338 }
1339 case 83: {
1340 double t = 2*y100 - 167;
1341 return 0.49243342227179841649e0 + (0.10037550043909497071e-1 + (0.11750334542845234952e-3 + (0.10544006716188967172e-5 + (0.73484461168242224872e-8 + (0.39383162326435752965e-10 + 0.15559069118222222222e-12 * t) * t) * t) * t) * t) * t;
1342 }
1343 case 84: {
1344 double t = 2*y100 - 169;
1345 return 0.51298708979209258326e0 + (0.10520454564612427224e-1 + (0.12400930037494996655e-3 + (0.11147886579371265246e-5 + (0.77517184550568711454e-8 + (0.41283980931872622611e-10 + 0.16122419680000000000e-12 * t) * t) * t) * t) * t) * t;
1346 }
1347 case 85: {
1348 double t = 2*y100 - 171;
1349 return 0.53453307979101369843e0 + (0.11030120618800726938e-1 + (0.13088741519572269581e-3 + (0.11784797595374515432e-5 + (0.81743383063044825400e-8 + (0.43252818449517081051e-10 + 0.16692592640000000000e-12 * t) * t) * t) * t) * t) * t;
1350 }
1351 case 86: {
1352 double t = 2*y100 - 173;
1353 return 0.55712643071169299478e0 + (0.11568077107929735233e-1 + (0.13815797838036651289e-3 + (0.12456314879260904558e-5 + (0.86169898078969313597e-8 + (0.45290446811539652525e-10 + 0.17268801084444444444e-12 * t) * t) * t) * t) * t) * t;
1354 }
1355 case 87: {
1356 double t = 2*y100 - 175;
1357 return 0.58082532122519320968e0 + (0.12135935999503877077e-1 + (0.14584223996665838559e-3 + (0.13164068573095710742e-5 + (0.90803643355106020163e-8 + (0.47397540713124619155e-10 + 0.17850211608888888889e-12 * t) * t) * t) * t) * t) * t;
1358 }
1359 case 88: {
1360 double t = 2*y100 - 177;
1361 return 0.60569124025293375554e0 + (0.12735396239525550361e-1 + (0.15396244472258863344e-3 + (0.13909744385382818253e-5 + (0.95651595032306228245e-8 + (0.49574672127669041550e-10 + 0.18435945564444444444e-12 * t) * t) * t) * t) * t) * t;
1362 }
1363 case 89: {
1364 double t = 2*y100 - 179;
1365 return 0.63178916494715716894e0 + (0.13368247798287030927e-1 + (0.16254186562762076141e-3 + (0.14695084048334056083e-5 + (0.10072078109604152350e-7 + (0.51822304995680707483e-10 + 0.19025081422222222222e-12 * t) * t) * t) * t) * t) * t;
1366 }
1367 case 90: {
1368 double t = 2*y100 - 181;
1369 return 0.65918774689725319200e0 + (0.14036375850601992063e-1 + (0.17160483760259706354e-3 + (0.15521885688723188371e-5 + (0.10601827031535280590e-7 + (0.54140790105837520499e-10 + 0.19616655146666666667e-12 * t) * t) * t) * t) * t) * t;
1370 }
1371 case 91: {
1372 double t = 2*y100 - 183;
1373 return 0.68795950683174433822e0 + (0.14741765091365869084e-1 + (0.18117679143520433835e-3 + (0.16392004108230585213e-5 + (0.11155116068018043001e-7 + (0.56530360194925690374e-10 + 0.20209663662222222222e-12 * t) * t) * t) * t) * t) * t;
1374 }
1375 case 92: {
1376 double t = 2*y100 - 185;
1377 return 0.71818103808729967036e0 + (0.15486504187117112279e-1 + (0.19128428784550923217e-3 + (0.17307350969359975848e-5 + (0.11732656736113607751e-7 + (0.58991125287563833603e-10 + 0.20803065333333333333e-12 * t) * t) * t) * t) * t) * t;
1378 }
1379 case 93: {
1380 double t = 2*y100 - 187;
1381 return 0.74993321911726254661e0 + (0.16272790364044783382e-1 + (0.20195505163377912645e-3 + (0.18269894883203346953e-5 + (0.12335161021630225535e-7 + (0.61523068312169087227e-10 + 0.21395783431111111111e-12 * t) * t) * t) * t) * t) * t;
1382 }
1383 case 94: {
1384 double t = 2*y100 - 189;
1385 return 0.78330143531283492729e0 + (0.17102934132652429240e-1 + (0.21321800585063327041e-3 + (0.19281661395543913713e-5 + (0.12963340087354341574e-7 + (0.64126040998066348872e-10 + 0.21986708942222222222e-12 * t) * t) * t) * t) * t) * t;
1386 }
1387 case 95: {
1388 double t = 2*y100 - 191;
1389 return 0.81837581041023811832e0 + (0.17979364149044223802e-1 + (0.22510330592753129006e-3 + (0.20344732868018175389e-5 + (0.13617902941839949718e-7 + (0.66799760083972474642e-10 + 0.22574701262222222222e-12 * t) * t) * t) * t) * t) * t;
1390 }
1391 case 96: {
1392 double t = 2*y100 - 193;
1393 return 0.85525144775685126237e0 + (0.18904632212547561026e-1 + (0.23764237370371255638e-3 + (0.21461248251306387979e-5 + (0.14299555071870523786e-7 + (0.69543803864694171934e-10 + 0.23158593688888888889e-12 * t) * t) * t) * t) * t) * t;
1394 }
1395 case 97: {
1396 double t = 2*y100 - 195;
1397 return 0.89402868170849933734e0 + (0.19881418399127202569e-1 + (0.25086793128395995798e-3 + (0.22633402747585233180e-5 + (0.15008997042116532283e-7 + (0.72357609075043941261e-10 + 0.23737194737777777778e-12 * t) * t) * t) * t) * t) * t;
1398 }
1399 case 98: {
1400 double t = 2*y100 - 197;
1401 return 0.93481333942870796363e0 + (0.20912536329780368893e-1 + (0.26481403465998477969e-3 + (0.23863447359754921676e-5 + (0.15746923065472184451e-7 + (0.75240468141720143653e-10 + 0.24309291271111111111e-12 * t) * t) * t) * t) * t) * t;
1402 }
1403 case 99: {
1404 double t = 2*y100 - 199;
1405 return 0.97771701335885035464e0 + (0.22000938572830479551e-1 + (0.27951610702682383001e-3 + (0.25153688325245314530e-5 + (0.16514019547822821453e-7 + (0.78191526829368231251e-10 + 0.24873652355555555556e-12 * t) * t) * t) * t) * t) * t;
1406 }
1407 }
1408 // we only get here if y = 1, i.e. |x| < 4*eps, in which case
1409 // erfcx is within 1e-15 of 1..
1410 return 1.0;
1411 }
1412
FADDEEVA_RE(erfcx)1413 double FADDEEVA_RE(erfcx)(double x)
1414 {
1415 if (x >= 0) {
1416 if (x > 50) { // continued-fraction expansion is faster
1417 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1418 if (x > 5e7) // 1-term expansion, important to avoid overflow
1419 return ispi / x;
1420 /* 5-term expansion (rely on compiler for CSE), simplified from:
1421 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */
1422 return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
1423 }
1424 return erfcx_y100(400/(4+x));
1425 }
1426 else
1427 return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x)
1428 : 2*exp(x*x) - erfcx_y100(400/(4-x)));
1429 }
1430
1431 /////////////////////////////////////////////////////////////////////////
1432 /* Compute a scaled Dawson integral
1433 FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
1434 equivalent to the imaginary part w(x) for real x.
1435
1436 Uses methods similar to the erfcx calculation above: continued fractions
1437 for large |x|, a lookup table of Chebyshev polynomials for smaller |x|,
1438 and finally a Taylor expansion for |x|<0.01.
1439
1440 Steven G. Johnson, October 2012. */
1441
1442 /* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).
1443
1444 Uses a look-up table of 100 different Chebyshev polynomials
1445 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
1446 with the help of Maple and a little shell script. This allows
1447 the Chebyshev polynomials to be of significantly lower degree (about 1/30)
1448 compared to fitting the whole [0,1] interval with a single polynomial. */
w_im_y100(double y100,double x)1449 static double w_im_y100(double y100, double x) {
1450 switch (static_cast<int> (y100)) {
1451 case 0: {
1452 double t = 2*y100 - 1;
1453 return 0.28351593328822191546e-2 + (0.28494783221378400759e-2 + (0.14427470563276734183e-4 + (0.10939723080231588129e-6 + (0.92474307943275042045e-9 + (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * t) * t) * t) * t) * t;
1454 }
1455 case 1: {
1456 double t = 2*y100 - 3;
1457 return 0.85927161243940350562e-2 + (0.29085312941641339862e-2 + (0.15106783707725582090e-4 + (0.11716709978531327367e-6 + (0.10197387816021040024e-8 + (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * t) * t) * t) * t) * t;
1458 }
1459 case 2: {
1460 double t = 2*y100 - 5;
1461 return 0.14471159831187703054e-1 + (0.29703978970263836210e-2 + (0.15835096760173030976e-4 + (0.12574803383199211596e-6 + (0.11278672159518415848e-8 + (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * t) * t) * t) * t) * t;
1462 }
1463 case 3: {
1464 double t = 2*y100 - 7;
1465 return 0.20476320420324610618e-1 + (0.30352843012898665856e-2 + (0.16617609387003727409e-4 + (0.13525429711163116103e-6 + (0.12515095552507169013e-8 + (0.13235687543603382345e-10 + 0.15326595042666666667e-12 * t) * t) * t) * t) * t) * t;
1466 }
1467 case 4: {
1468 double t = 2*y100 - 9;
1469 return 0.26614461952489004566e-1 + (0.31034189276234947088e-2 + (0.17460268109986214274e-4 + (0.14582130824485709573e-6 + (0.13935959083809746345e-8 + (0.15249438072998932900e-10 + 0.18344741882133333333e-12 * t) * t) * t) * t) * t) * t;
1470 }
1471 case 5: {
1472 double t = 2*y100 - 11;
1473 return 0.32892330248093586215e-1 + (0.31750557067975068584e-2 + (0.18369907582308672632e-4 + (0.15761063702089457882e-6 + (0.15577638230480894382e-8 + (0.17663868462699097951e-10 + (0.22126732680711111111e-12 + 0.30273474177737853668e-14 * t) * t) * t) * t) * t) * t) * t;
1474 }
1475 case 6: {
1476 double t = 2*y100 - 13;
1477 return 0.39317207681134336024e-1 + (0.32504779701937539333e-2 + (0.19354426046513400534e-4 + (0.17081646971321290539e-6 + (0.17485733959327106250e-8 + (0.20593687304921961410e-10 + (0.26917401949155555556e-12 + 0.38562123837725712270e-14 * t) * t) * t) * t) * t) * t) * t;
1478 }
1479 case 7: {
1480 double t = 2*y100 - 15;
1481 return 0.45896976511367738235e-1 + (0.33300031273110976165e-2 + (0.20423005398039037313e-4 + (0.18567412470376467303e-6 + (0.19718038363586588213e-8 + (0.24175006536781219807e-10 + (0.33059982791466666666e-12 + 0.49756574284439426165e-14 * t) * t) * t) * t) * t) * t) * t;
1482 }
1483 case 8: {
1484 double t = 2*y100 - 17;
1485 return 0.52640192524848962855e-1 + (0.34139883358846720806e-2 + (0.21586390240603337337e-4 + (0.20247136501568904646e-6 + (0.22348696948197102935e-8 + (0.28597516301950162548e-10 + (0.41045502119111111110e-12 + 0.65151614515238361946e-14 * t) * t) * t) * t) * t) * t) * t;
1486 }
1487 case 9: {
1488 double t = 2*y100 - 19;
1489 return 0.59556171228656770456e-1 + (0.35028374386648914444e-2 + (0.22857246150998562824e-4 + (0.22156372146525190679e-6 + (0.25474171590893813583e-8 + (0.34122390890697400584e-10 + (0.51593189879111111110e-12 + 0.86775076853908006938e-14 * t) * t) * t) * t) * t) * t) * t;
1490 }
1491 case 10: {
1492 double t = 2*y100 - 21;
1493 return 0.66655089485108212551e-1 + (0.35970095381271285568e-2 + (0.24250626164318672928e-4 + (0.24339561521785040536e-6 + (0.29221990406518411415e-8 + (0.41117013527967776467e-10 + (0.65786450716444444445e-12 + 0.11791885745450623331e-13 * t) * t) * t) * t) * t) * t) * t;
1494 }
1495 case 11: {
1496 double t = 2*y100 - 23;
1497 return 0.73948106345519174661e-1 + (0.36970297216569341748e-2 + (0.25784588137312868792e-4 + (0.26853012002366752770e-6 + (0.33763958861206729592e-8 + (0.50111549981376976397e-10 + (0.85313857496888888890e-12 + 0.16417079927706899860e-13 * t) * t) * t) * t) * t) * t) * t;
1498 }
1499 case 12: {
1500 double t = 2*y100 - 25;
1501 return 0.81447508065002963203e-1 + (0.38035026606492705117e-2 + (0.27481027572231851896e-4 + (0.29769200731832331364e-6 + (0.39336816287457655076e-8 + (0.61895471132038157624e-10 + (0.11292303213511111111e-11 + 0.23558532213703884304e-13 * t) * t) * t) * t) * t) * t) * t;
1502 }
1503 case 13: {
1504 double t = 2*y100 - 27;
1505 return 0.89166884027582716628e-1 + (0.39171301322438946014e-2 + (0.29366827260422311668e-4 + (0.33183204390350724895e-6 + (0.46276006281647330524e-8 + (0.77692631378169813324e-10 + (0.15335153258844444444e-11 + 0.35183103415916026911e-13 * t) * t) * t) * t) * t) * t) * t;
1506 }
1507 case 14: {
1508 double t = 2*y100 - 29;
1509 return 0.97121342888032322019e-1 + (0.40387340353207909514e-2 + (0.31475490395950776930e-4 + (0.37222714227125135042e-6 + (0.55074373178613809996e-8 + (0.99509175283990337944e-10 + (0.21552645758222222222e-11 + 0.55728651431872687605e-13 * t) * t) * t) * t) * t) * t) * t;
1510 }
1511 case 15: {
1512 double t = 2*y100 - 31;
1513 return 0.10532778218603311137e0 + (0.41692873614065380607e-2 + (0.33849549774889456984e-4 + (0.42064596193692630143e-6 + (0.66494579697622432987e-8 + (0.13094103581931802337e-9 + (0.31896187409777777778e-11 + 0.97271974184476560742e-13 * t) * t) * t) * t) * t) * t) * t;
1514 }
1515 case 16: {
1516 double t = 2*y100 - 33;
1517 return 0.11380523107427108222e0 + (0.43099572287871821013e-2 + (0.36544324341565929930e-4 + (0.47965044028581857764e-6 + (0.81819034238463698796e-8 + (0.17934133239549647357e-9 + (0.50956666166186293627e-11 + (0.18850487318190638010e-12 + 0.79697813173519853340e-14 * t) * t) * t) * t) * t) * t) * t) * t;
1518 }
1519 case 17: {
1520 double t = 2*y100 - 35;
1521 return 0.12257529703447467345e0 + (0.44621675710026986366e-2 + (0.39634304721292440285e-4 + (0.55321553769873381819e-6 + (0.10343619428848520870e-7 + (0.26033830170470368088e-9 + (0.87743837749108025357e-11 + (0.34427092430230063401e-12 + 0.10205506615709843189e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1522 }
1523 case 18: {
1524 double t = 2*y100 - 37;
1525 return 0.13166276955656699478e0 + (0.46276970481783001803e-2 + (0.43225026380496399310e-4 + (0.64799164020016902656e-6 + (0.13580082794704641782e-7 + (0.39839800853954313927e-9 + (0.14431142411840000000e-10 + 0.42193457308830027541e-12 * t) * t) * t) * t) * t) * t) * t;
1526 }
1527 case 19: {
1528 double t = 2*y100 - 39;
1529 return 0.14109647869803356475e0 + (0.48088424418545347758e-2 + (0.47474504753352150205e-4 + (0.77509866468724360352e-6 + (0.18536851570794291724e-7 + (0.60146623257887570439e-9 + (0.18533978397305276318e-10 + (0.41033845938901048380e-13 - 0.46160680279304825485e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1530 }
1531 case 20: {
1532 double t = 2*y100 - 41;
1533 return 0.15091057940548936603e0 + (0.50086864672004685703e-2 + (0.52622482832192230762e-4 + (0.95034664722040355212e-6 + (0.25614261331144718769e-7 + (0.80183196716888606252e-9 + (0.12282524750534352272e-10 + (-0.10531774117332273617e-11 - 0.86157181395039646412e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1534 }
1535 case 21: {
1536 double t = 2*y100 - 43;
1537 return 0.16114648116017010770e0 + (0.52314661581655369795e-2 + (0.59005534545908331315e-4 + (0.11885518333915387760e-5 + (0.33975801443239949256e-7 + (0.82111547144080388610e-9 + (-0.12357674017312854138e-10 + (-0.24355112256914479176e-11 - 0.75155506863572930844e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1538 }
1539 case 22: {
1540 double t = 2*y100 - 45;
1541 return 0.17185551279680451144e0 + (0.54829002967599420860e-2 + (0.67013226658738082118e-4 + (0.14897400671425088807e-5 + (0.40690283917126153701e-7 + (0.44060872913473778318e-9 + (-0.52641873433280000000e-10 - 0.30940587864543343124e-11 * t) * t) * t) * t) * t) * t) * t;
1542 }
1543 case 23: {
1544 double t = 2*y100 - 47;
1545 return 0.18310194559815257381e0 + (0.57701559375966953174e-2 + (0.76948789401735193483e-4 + (0.18227569842290822512e-5 + (0.41092208344387212276e-7 + (-0.44009499965694442143e-9 + (-0.92195414685628803451e-10 + (-0.22657389705721753299e-11 + 0.10004784908106839254e-12 * t) * t) * t) * t) * t) * t) * t) * t;
1546 }
1547 case 24: {
1548 double t = 2*y100 - 49;
1549 return 0.19496527191546630345e0 + (0.61010853144364724856e-2 + (0.88812881056342004864e-4 + (0.21180686746360261031e-5 + (0.30652145555130049203e-7 + (-0.16841328574105890409e-8 + (-0.11008129460612823934e-9 + (-0.12180794204544515779e-12 + 0.15703325634590334097e-12 * t) * t) * t) * t) * t) * t) * t) * t;
1550 }
1551 case 25: {
1552 double t = 2*y100 - 51;
1553 return 0.20754006813966575720e0 + (0.64825787724922073908e-2 + (0.10209599627522311893e-3 + (0.22785233392557600468e-5 + (0.73495224449907568402e-8 + (-0.29442705974150112783e-8 + (-0.94082603434315016546e-10 + (0.23609990400179321267e-11 + 0.14141908654269023788e-12 * t) * t) * t) * t) * t) * t) * t) * t;
1554 }
1555 case 26: {
1556 double t = 2*y100 - 53;
1557 return 0.22093185554845172146e0 + (0.69182878150187964499e-2 + (0.11568723331156335712e-3 + (0.22060577946323627739e-5 + (-0.26929730679360840096e-7 + (-0.38176506152362058013e-8 + (-0.47399503861054459243e-10 + (0.40953700187172127264e-11 + 0.69157730376118511127e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1558 }
1559 case 27: {
1560 double t = 2*y100 - 55;
1561 return 0.23524827304057813918e0 + (0.74063350762008734520e-2 + (0.12796333874615790348e-3 + (0.18327267316171054273e-5 + (-0.66742910737957100098e-7 + (-0.40204740975496797870e-8 + (0.14515984139495745330e-10 + (0.44921608954536047975e-11 - 0.18583341338983776219e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1562 }
1563 case 28: {
1564 double t = 2*y100 - 57;
1565 return 0.25058626331812744775e0 + (0.79377285151602061328e-2 + (0.13704268650417478346e-3 + (0.11427511739544695861e-5 + (-0.10485442447768377485e-6 + (-0.34850364756499369763e-8 + (0.72656453829502179208e-10 + (0.36195460197779299406e-11 - 0.84882136022200714710e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1566 }
1567 case 29: {
1568 double t = 2*y100 - 59;
1569 return 0.26701724900280689785e0 + (0.84959936119625864274e-2 + (0.14112359443938883232e-3 + (0.17800427288596909634e-6 + (-0.13443492107643109071e-6 + (-0.23512456315677680293e-8 + (0.11245846264695936769e-9 + (0.19850501334649565404e-11 - 0.11284666134635050832e-12 * t) * t) * t) * t) * t) * t) * t) * t;
1570 }
1571 case 30: {
1572 double t = 2*y100 - 61;
1573 return 0.28457293586253654144e0 + (0.90581563892650431899e-2 + (0.13880520331140646738e-3 + (-0.97262302362522896157e-6 + (-0.15077100040254187366e-6 + (-0.88574317464577116689e-9 + (0.12760311125637474581e-9 + (0.20155151018282695055e-12 - 0.10514169375181734921e-12 * t) * t) * t) * t) * t) * t) * t) * t;
1574 }
1575 case 31: {
1576 double t = 2*y100 - 63;
1577 return 0.30323425595617385705e0 + (0.95968346790597422934e-2 + (0.12931067776725883939e-3 + (-0.21938741702795543986e-5 + (-0.15202888584907373963e-6 + (0.61788350541116331411e-9 + (0.11957835742791248256e-9 + (-0.12598179834007710908e-11 - 0.75151817129574614194e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1578 }
1579 case 32: {
1580 double t = 2*y100 - 65;
1581 return 0.32292521181517384379e0 + (0.10082957727001199408e-1 + (0.11257589426154962226e-3 + (-0.33670890319327881129e-5 + (-0.13910529040004008158e-6 + (0.19170714373047512945e-8 + (0.94840222377720494290e-10 + (-0.21650018351795353201e-11 - 0.37875211678024922689e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1582 }
1583 case 33: {
1584 double t = 2*y100 - 67;
1585 return 0.34351233557911753862e0 + (0.10488575435572745309e-1 + (0.89209444197248726614e-4 + (-0.43893459576483345364e-5 + (-0.11488595830450424419e-6 + (0.28599494117122464806e-8 + (0.61537542799857777779e-10 - 0.24935749227658002212e-11 * t) * t) * t) * t) * t) * t) * t;
1586 }
1587 case 34: {
1588 double t = 2*y100 - 69;
1589 return 0.36480946642143669093e0 + (0.10789304203431861366e-1 + (0.60357993745283076834e-4 + (-0.51855862174130669389e-5 + (-0.83291664087289801313e-7 + (0.33898011178582671546e-8 + (0.27082948188277716482e-10 + (-0.23603379397408694974e-11 + 0.19328087692252869842e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1590 }
1591 case 35: {
1592 double t = 2*y100 - 71;
1593 return 0.38658679935694939199e0 + (0.10966119158288804999e-1 + (0.27521612041849561426e-4 + (-0.57132774537670953638e-5 + (-0.48404772799207914899e-7 + (0.35268354132474570493e-8 + (-0.32383477652514618094e-11 + (-0.19334202915190442501e-11 + 0.32333189861286460270e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1594 }
1595 case 36: {
1596 double t = 2*y100 - 73;
1597 return 0.40858275583808707870e0 + (0.11006378016848466550e-1 + (-0.76396376685213286033e-5 + (-0.59609835484245791439e-5 + (-0.13834610033859313213e-7 + (0.33406952974861448790e-8 + (-0.26474915974296612559e-10 + (-0.13750229270354351983e-11 + 0.36169366979417390637e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1598 }
1599 case 37: {
1600 double t = 2*y100 - 75;
1601 return 0.43051714914006682977e0 + (0.10904106549500816155e-1 + (-0.43477527256787216909e-4 + (-0.59429739547798343948e-5 + (0.17639200194091885949e-7 + (0.29235991689639918688e-8 + (-0.41718791216277812879e-10 + (-0.81023337739508049606e-12 + 0.33618915934461994428e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1602 }
1603 case 38: {
1604 double t = 2*y100 - 77;
1605 return 0.45210428135559607406e0 + (0.10659670756384400554e-1 + (-0.78488639913256978087e-4 + (-0.56919860886214735936e-5 + (0.44181850467477733407e-7 + (0.23694306174312688151e-8 + (-0.49492621596685443247e-10 + (-0.31827275712126287222e-12 + 0.27494438742721623654e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1606 }
1607 case 39: {
1608 double t = 2*y100 - 79;
1609 return 0.47306491195005224077e0 + (0.10279006119745977570e-1 + (-0.11140268171830478306e-3 + (-0.52518035247451432069e-5 + (0.64846898158889479518e-7 + (0.17603624837787337662e-8 + (-0.51129481592926104316e-10 + (0.62674584974141049511e-13 + 0.20055478560829935356e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1610 }
1611 case 40: {
1612 double t = 2*y100 - 81;
1613 return 0.49313638965719857647e0 + (0.97725799114772017662e-2 + (-0.14122854267291533334e-3 + (-0.46707252568834951907e-5 + (0.79421347979319449524e-7 + (0.11603027184324708643e-8 + (-0.48269605844397175946e-10 + (0.32477251431748571219e-12 + 0.12831052634143527985e-13 * t) * t) * t) * t) * t) * t) * t) * t;
1614 }
1615 case 41: {
1616 double t = 2*y100 - 83;
1617 return 0.51208057433416004042e0 + (0.91542422354009224951e-2 + (-0.16726530230228647275e-3 + (-0.39964621752527649409e-5 + (0.88232252903213171454e-7 + (0.61343113364949928501e-9 + (-0.42516755603130443051e-10 + (0.47910437172240209262e-12 + 0.66784341874437478953e-14 * t) * t) * t) * t) * t) * t) * t) * t;
1618 }
1619 case 42: {
1620 double t = 2*y100 - 85;
1621 return 0.52968945458607484524e0 + (0.84400880445116786088e-2 + (-0.18908729783854258774e-3 + (-0.32725905467782951931e-5 + (0.91956190588652090659e-7 + (0.14593989152420122909e-9 + (-0.35239490687644444445e-10 + 0.54613829888448694898e-12 * t) * t) * t) * t) * t) * t) * t;
1622 }
1623 case 43: {
1624 double t = 2*y100 - 87;
1625 return 0.54578857454330070965e0 + (0.76474155195880295311e-2 + (-0.20651230590808213884e-3 + (-0.25364339140543131706e-5 + (0.91455367999510681979e-7 + (-0.23061359005297528898e-9 + (-0.27512928625244444444e-10 + 0.54895806008493285579e-12 * t) * t) * t) * t) * t) * t) * t;
1626 }
1627 case 44: {
1628 double t = 2*y100 - 89;
1629 return 0.56023851910298493910e0 + (0.67938321739997196804e-2 + (-0.21956066613331411760e-3 + (-0.18181127670443266395e-5 + (0.87650335075416845987e-7 + (-0.51548062050366615977e-9 + (-0.20068462174044444444e-10 + 0.50912654909758187264e-12 * t) * t) * t) * t) * t) * t) * t;
1630 }
1631 case 45: {
1632 double t = 2*y100 - 91;
1633 return 0.57293478057455721150e0 + (0.58965321010394044087e-2 + (-0.22841145229276575597e-3 + (-0.11404605562013443659e-5 + (0.81430290992322326296e-7 + (-0.71512447242755357629e-9 + (-0.13372664928000000000e-10 + 0.44461498336689298148e-12 * t) * t) * t) * t) * t) * t) * t;
1634 }
1635 case 46: {
1636 double t = 2*y100 - 93;
1637 return 0.58380635448407827360e0 + (0.49717469530842831182e-2 + (-0.23336001540009645365e-3 + (-0.51952064448608850822e-6 + (0.73596577815411080511e-7 + (-0.84020916763091566035e-9 + (-0.76700972702222222221e-11 + 0.36914462807972467044e-12 * t) * t) * t) * t) * t) * t) * t;
1638 }
1639 case 47: {
1640 double t = 2*y100 - 95;
1641 return 0.59281340237769489597e0 + (0.40343592069379730568e-2 + (-0.23477963738658326185e-3 + (0.34615944987790224234e-7 + (0.64832803248395814574e-7 + (-0.90329163587627007971e-9 + (-0.30421940400000000000e-11 + 0.29237386653743536669e-12 * t) * t) * t) * t) * t) * t) * t;
1642 }
1643 case 48: {
1644 double t = 2*y100 - 97;
1645 return 0.59994428743114271918e0 + (0.30976579788271744329e-2 + (-0.23308875765700082835e-3 + (0.51681681023846925160e-6 + (0.55694594264948268169e-7 + (-0.91719117313243464652e-9 + (0.53982743680000000000e-12 + 0.22050829296187771142e-12 * t) * t) * t) * t) * t) * t) * t;
1646 }
1647 case 49: {
1648 double t = 2*y100 - 99;
1649 return 0.60521224471819875444e0 + (0.21732138012345456060e-2 + (-0.22872428969625997456e-3 + (0.92588959922653404233e-6 + (0.46612665806531930684e-7 + (-0.89393722514414153351e-9 + (0.31718550353777777778e-11 + 0.15705458816080549117e-12 * t) * t) * t) * t) * t) * t) * t;
1650 }
1651 case 50: {
1652 double t = 2*y100 - 101;
1653 return 0.60865189969791123620e0 + (0.12708480848877451719e-2 + (-0.22212090111534847166e-3 + (0.12636236031532793467e-5 + (0.37904037100232937574e-7 + (-0.84417089968101223519e-9 + (0.49843180828444444445e-11 + 0.10355439441049048273e-12 * t) * t) * t) * t) * t) * t) * t;
1654 }
1655 case 51: {
1656 double t = 2*y100 - 103;
1657 return 0.61031580103499200191e0 + (0.39867436055861038223e-3 + (-0.21369573439579869291e-3 + (0.15339402129026183670e-5 + (0.29787479206646594442e-7 + (-0.77687792914228632974e-9 + (0.61192452741333333334e-11 + 0.60216691829459295780e-13 * t) * t) * t) * t) * t) * t) * t;
1658 }
1659 case 52: {
1660 double t = 2*y100 - 105;
1661 return 0.61027109047879835868e0 + (-0.43680904508059878254e-3 + (-0.20383783788303894442e-3 + (0.17421743090883439959e-5 + (0.22400425572175715576e-7 + (-0.69934719320045128997e-9 + (0.67152759655111111110e-11 + 0.26419960042578359995e-13 * t) * t) * t) * t) * t) * t) * t;
1662 }
1663 case 53: {
1664 double t = 2*y100 - 107;
1665 return 0.60859639489217430521e0 + (-0.12305921390962936873e-2 + (-0.19290150253894682629e-3 + (0.18944904654478310128e-5 + (0.15815530398618149110e-7 + (-0.61726850580964876070e-9 + 0.68987888999111111110e-11 * t) * t) * t) * t) * t) * t;
1666 }
1667 case 54: {
1668 double t = 2*y100 - 109;
1669 return 0.60537899426486075181e0 + (-0.19790062241395705751e-2 + (-0.18120271393047062253e-3 + (0.19974264162313241405e-5 + (0.10055795094298172492e-7 + (-0.53491997919318263593e-9 + (0.67794550295111111110e-11 - 0.17059208095741511603e-13 * t) * t) * t) * t) * t) * t) * t;
1670 }
1671 case 55: {
1672 double t = 2*y100 - 111;
1673 return 0.60071229457904110537e0 + (-0.26795676776166354354e-2 + (-0.16901799553627508781e-3 + (0.20575498324332621581e-5 + (0.51077165074461745053e-8 + (-0.45536079828057221858e-9 + (0.64488005516444444445e-11 - 0.29311677573152766338e-13 * t) * t) * t) * t) * t) * t) * t;
1674 }
1675 case 56: {
1676 double t = 2*y100 - 113;
1677 return 0.59469361520112714738e0 + (-0.33308208190600993470e-2 + (-0.15658501295912405679e-3 + (0.20812116912895417272e-5 + (0.93227468760614182021e-9 + (-0.38066673740116080415e-9 + (0.59806790359111111110e-11 - 0.36887077278950440597e-13 * t) * t) * t) * t) * t) * t) * t;
1678 }
1679 case 57: {
1680 double t = 2*y100 - 115;
1681 return 0.58742228631775388268e0 + (-0.39321858196059227251e-2 + (-0.14410441141450122535e-3 + (0.20743790018404020716e-5 + (-0.25261903811221913762e-8 + (-0.31212416519526924318e-9 + (0.54328422462222222221e-11 - 0.40864152484979815972e-13 * t) * t) * t) * t) * t) * t) * t;
1682 }
1683 case 58: {
1684 double t = 2*y100 - 117;
1685 return 0.57899804200033018447e0 + (-0.44838157005618913447e-2 + (-0.13174245966501437965e-3 + (0.20425306888294362674e-5 + (-0.53330296023875447782e-8 + (-0.25041289435539821014e-9 + (0.48490437205333333334e-11 - 0.42162206939169045177e-13 * t) * t) * t) * t) * t) * t) * t;
1686 }
1687 case 59: {
1688 double t = 2*y100 - 119;
1689 return 0.56951968796931245974e0 + (-0.49864649488074868952e-2 + (-0.11963416583477567125e-3 + (0.19906021780991036425e-5 + (-0.75580140299436494248e-8 + (-0.19576060961919820491e-9 + (0.42613011928888888890e-11 - 0.41539443304115604377e-13 * t) * t) * t) * t) * t) * t) * t;
1690 }
1691 case 60: {
1692 double t = 2*y100 - 121;
1693 return 0.55908401930063918964e0 + (-0.54413711036826877753e-2 + (-0.10788661102511914628e-3 + (0.19229663322982839331e-5 + (-0.92714731195118129616e-8 + (-0.14807038677197394186e-9 + (0.36920870298666666666e-11 - 0.39603726688419162617e-13 * t) * t) * t) * t) * t) * t) * t;
1694 }
1695 case 61: {
1696 double t = 2*y100 - 123;
1697 return 0.54778496152925675315e0 + (-0.58501497933213396670e-2 + (-0.96582314317855227421e-4 + (0.18434405235069270228e-5 + (-0.10541580254317078711e-7 + (-0.10702303407788943498e-9 + (0.31563175582222222222e-11 - 0.36829748079110481422e-13 * t) * t) * t) * t) * t) * t) * t;
1698 }
1699 case 62: {
1700 double t = 2*y100 - 125;
1701 return 0.53571290831682823999e0 + (-0.62147030670760791791e-2 + (-0.85782497917111760790e-4 + (0.17553116363443470478e-5 + (-0.11432547349815541084e-7 + (-0.72157091369041330520e-10 + (0.26630811607111111111e-11 - 0.33578660425893164084e-13 * t) * t) * t) * t) * t) * t) * t;
1702 }
1703 case 63: {
1704 double t = 2*y100 - 127;
1705 return 0.52295422962048434978e0 + (-0.65371404367776320720e-2 + (-0.75530164941473343780e-4 + (0.16613725797181276790e-5 + (-0.12003521296598910761e-7 + (-0.42929753689181106171e-10 + (0.22170894940444444444e-11 - 0.30117697501065110505e-13 * t) * t) * t) * t) * t) * t) * t;
1706 }
1707 case 64: {
1708 double t = 2*y100 - 129;
1709 return 0.50959092577577886140e0 + (-0.68197117603118591766e-2 + (-0.65852936198953623307e-4 + (0.15639654113906716939e-5 + (-0.12308007991056524902e-7 + (-0.18761997536910939570e-10 + (0.18198628922666666667e-11 - 0.26638355362285200932e-13 * t) * t) * t) * t) * t) * t) * t;
1710 }
1711 case 65: {
1712 double t = 2*y100 - 131;
1713 return 0.49570040481823167970e0 + (-0.70647509397614398066e-2 + (-0.56765617728962588218e-4 + (0.14650274449141448497e-5 + (-0.12393681471984051132e-7 + (0.92904351801168955424e-12 + (0.14706755960177777778e-11 - 0.23272455351266325318e-13 * t) * t) * t) * t) * t) * t) * t;
1714 }
1715 case 66: {
1716 double t = 2*y100 - 133;
1717 return 0.48135536250935238066e0 + (-0.72746293327402359783e-2 + (-0.48272489495730030780e-4 + (0.13661377309113939689e-5 + (-0.12302464447599382189e-7 + (0.16707760028737074907e-10 + (0.11672928324444444444e-11 - 0.20105801424709924499e-13 * t) * t) * t) * t) * t) * t) * t;
1718 }
1719 case 67: {
1720 double t = 2*y100 - 135;
1721 return 0.46662374675511439448e0 + (-0.74517177649528487002e-2 + (-0.40369318744279128718e-4 + (0.12685621118898535407e-5 + (-0.12070791463315156250e-7 + (0.29105507892605823871e-10 + (0.90653314645333333334e-12 - 0.17189503312102982646e-13 * t) * t) * t) * t) * t) * t) * t;
1722 }
1723 case 68: {
1724 double t = 2*y100 - 137;
1725 return 0.45156879030168268778e0 + (-0.75983560650033817497e-2 + (-0.33045110380705139759e-4 + (0.11732956732035040896e-5 + (-0.11729986947158201869e-7 + (0.38611905704166441308e-10 + (0.68468768305777777779e-12 - 0.14549134330396754575e-13 * t) * t) * t) * t) * t) * t) * t;
1726 }
1727 case 69: {
1728 double t = 2*y100 - 139;
1729 return 0.43624909769330896904e0 + (-0.77168291040309554679e-2 + (-0.26283612321339907756e-4 + (0.10811018836893550820e-5 + (-0.11306707563739851552e-7 + (0.45670446788529607380e-10 + (0.49782492549333333334e-12 - 0.12191983967561779442e-13 * t) * t) * t) * t) * t) * t) * t;
1730 }
1731 case 70: {
1732 double t = 2*y100 - 141;
1733 return 0.42071877443548481181e0 + (-0.78093484015052730097e-2 + (-0.20064596897224934705e-4 + (0.99254806680671890766e-6 + (-0.10823412088884741451e-7 + (0.50677203326904716247e-10 + (0.34200547594666666666e-12 - 0.10112698698356194618e-13 * t) * t) * t) * t) * t) * t) * t;
1734 }
1735 case 71: {
1736 double t = 2*y100 - 143;
1737 return 0.40502758809710844280e0 + (-0.78780384460872937555e-2 + (-0.14364940764532853112e-4 + (0.90803709228265217384e-6 + (-0.10298832847014466907e-7 + (0.53981671221969478551e-10 + (0.21342751381333333333e-12 - 0.82975901848387729274e-14 * t) * t) * t) * t) * t) * t) * t;
1738 }
1739 case 72: {
1740 double t = 2*y100 - 145;
1741 return 0.38922115269731446690e0 + (-0.79249269708242064120e-2 + (-0.91595258799106970453e-5 + (0.82783535102217576495e-6 + (-0.97484311059617744437e-8 + (0.55889029041660225629e-10 + (0.10851981336888888889e-12 - 0.67278553237853459757e-14 * t) * t) * t) * t) * t) * t) * t;
1742 }
1743 case 73: {
1744 double t = 2*y100 - 147;
1745 return 0.37334112915460307335e0 + (-0.79519385109223148791e-2 + (-0.44219833548840469752e-5 + (0.75209719038240314732e-6 + (-0.91848251458553190451e-8 + (0.56663266668051433844e-10 + (0.23995894257777777778e-13 - 0.53819475285389344313e-14 * t) * t) * t) * t) * t) * t) * t;
1746 }
1747 case 74: {
1748 double t = 2*y100 - 149;
1749 return 0.35742543583374223085e0 + (-0.79608906571527956177e-2 + (-0.12530071050975781198e-6 + (0.68088605744900552505e-6 + (-0.86181844090844164075e-8 + (0.56530784203816176153e-10 + (-0.43120012248888888890e-13 - 0.42372603392496813810e-14 * t) * t) * t) * t) * t) * t) * t;
1750 }
1751 case 75: {
1752 double t = 2*y100 - 151;
1753 return 0.34150846431979618536e0 + (-0.79534924968773806029e-2 + (0.37576885610891515813e-5 + (0.61419263633090524326e-6 + (-0.80565865409945960125e-8 + (0.55684175248749269411e-10 + (-0.95486860764444444445e-13 - 0.32712946432984510595e-14 * t) * t) * t) * t) * t) * t) * t;
1754 }
1755 case 76: {
1756 double t = 2*y100 - 153;
1757 return 0.32562129649136346824e0 + (-0.79313448067948884309e-2 + (0.72539159933545300034e-5 + (0.55195028297415503083e-6 + (-0.75063365335570475258e-8 + (0.54281686749699595941e-10 - 0.13545424295111111111e-12 * t) * t) * t) * t) * t) * t;
1758 }
1759 case 77: {
1760 double t = 2*y100 - 155;
1761 return 0.30979191977078391864e0 + (-0.78959416264207333695e-2 + (0.10389774377677210794e-4 + (0.49404804463196316464e-6 + (-0.69722488229411164685e-8 + (0.52469254655951393842e-10 - 0.16507860650666666667e-12 * t) * t) * t) * t) * t) * t;
1762 }
1763 case 78: {
1764 double t = 2*y100 - 157;
1765 return 0.29404543811214459904e0 + (-0.78486728990364155356e-2 + (0.13190885683106990459e-4 + (0.44034158861387909694e-6 + (-0.64578942561562616481e-8 + (0.50354306498006928984e-10 - 0.18614473550222222222e-12 * t) * t) * t) * t) * t) * t;
1766 }
1767 case 79: {
1768 double t = 2*y100 - 159;
1769 return 0.27840427686253660515e0 + (-0.77908279176252742013e-2 + (0.15681928798708548349e-4 + (0.39066226205099807573e-6 + (-0.59658144820660420814e-8 + (0.48030086420373141763e-10 - 0.20018995173333333333e-12 * t) * t) * t) * t) * t) * t;
1770 }
1771 case 80: {
1772 double t = 2*y100 - 161;
1773 return 0.26288838011163800908e0 + (-0.77235993576119469018e-2 + (0.17886516796198660969e-4 + (0.34482457073472497720e-6 + (-0.54977066551955420066e-8 + (0.45572749379147269213e-10 - 0.20852924954666666667e-12 * t) * t) * t) * t) * t) * t;
1774 }
1775 case 81: {
1776 double t = 2*y100 - 163;
1777 return 0.24751539954181029717e0 + (-0.76480877165290370975e-2 + (0.19827114835033977049e-4 + (0.30263228619976332110e-6 + (-0.50545814570120129947e-8 + (0.43043879374212005966e-10 - 0.21228012028444444444e-12 * t) * t) * t) * t) * t) * t;
1778 }
1779 case 82: {
1780 double t = 2*y100 - 165;
1781 return 0.23230087411688914593e0 + (-0.75653060136384041587e-2 + (0.21524991113020016415e-4 + (0.26388338542539382413e-6 + (-0.46368974069671446622e-8 + (0.40492715758206515307e-10 - 0.21238627815111111111e-12 * t) * t) * t) * t) * t) * t;
1782 }
1783 case 83: {
1784 double t = 2*y100 - 167;
1785 return 0.21725840021297341931e0 + (-0.74761846305979730439e-2 + (0.23000194404129495243e-4 + (0.22837400135642906796e-6 + (-0.42446743058417541277e-8 + (0.37958104071765923728e-10 - 0.20963978568888888889e-12 * t) * t) * t) * t) * t) * t;
1786 }
1787 case 84: {
1788 double t = 2*y100 - 169;
1789 return 0.20239979200788191491e0 + (-0.73815761980493466516e-2 + (0.24271552727631854013e-4 + (0.19590154043390012843e-6 + (-0.38775884642456551753e-8 + (0.35470192372162901168e-10 - 0.20470131678222222222e-12 * t) * t) * t) * t) * t) * t;
1790 }
1791 case 85: {
1792 double t = 2*y100 - 171;
1793 return 0.18773523211558098962e0 + (-0.72822604530339834448e-2 + (0.25356688567841293697e-4 + (0.16626710297744290016e-6 + (-0.35350521468015310830e-8 + (0.33051896213898864306e-10 - 0.19811844544000000000e-12 * t) * t) * t) * t) * t) * t;
1794 }
1795 case 86: {
1796 double t = 2*y100 - 173;
1797 return 0.17327341258479649442e0 + (-0.71789490089142761950e-2 + (0.26272046822383820476e-4 + (0.13927732375657362345e-6 + (-0.32162794266956859603e-8 + (0.30720156036105652035e-10 - 0.19034196304000000000e-12 * t) * t) * t) * t) * t) * t;
1798 }
1799 case 87: {
1800 double t = 2*y100 - 175;
1801 return 0.15902166648328672043e0 + (-0.70722899934245504034e-2 + (0.27032932310132226025e-4 + (0.11474573347816568279e-6 + (-0.29203404091754665063e-8 + (0.28487010262547971859e-10 - 0.18174029063111111111e-12 * t) * t) * t) * t) * t) * t;
1802 }
1803 case 88: {
1804 double t = 2*y100 - 177;
1805 return 0.14498609036610283865e0 + (-0.69628725220045029273e-2 + (0.27653554229160596221e-4 + (0.92493727167393036470e-7 + (-0.26462055548683583849e-8 + (0.26360506250989943739e-10 - 0.17261211260444444444e-12 * t) * t) * t) * t) * t) * t;
1806 }
1807 case 89: {
1808 double t = 2*y100 - 179;
1809 return 0.13117165798208050667e0 + (-0.68512309830281084723e-2 + (0.28147075431133863774e-4 + (0.72351212437979583441e-7 + (-0.23927816200314358570e-8 + (0.24345469651209833155e-10 - 0.16319736960000000000e-12 * t) * t) * t) * t) * t) * t;
1810 }
1811 case 90: {
1812 double t = 2*y100 - 181;
1813 return 0.11758232561160626306e0 + (-0.67378491192463392927e-2 + (0.28525664781722907847e-4 + (0.54156999310046790024e-7 + (-0.21589405340123827823e-8 + (0.22444150951727334619e-10 - 0.15368675584000000000e-12 * t) * t) * t) * t) * t) * t;
1814 }
1815 case 91: {
1816 double t = 2*y100 - 183;
1817 return 0.10422112945361673560e0 + (-0.66231638959845581564e-2 + (0.28800551216363918088e-4 + (0.37758983397952149613e-7 + (-0.19435423557038933431e-8 + (0.20656766125421362458e-10 - 0.14422990012444444444e-12 * t) * t) * t) * t) * t) * t;
1818 }
1819 case 92: {
1820 double t = 2*y100 - 185;
1821 return 0.91090275493541084785e-1 + (-0.65075691516115160062e-2 + (0.28982078385527224867e-4 + (0.23014165807643012781e-7 + (-0.17454532910249875958e-8 + (0.18981946442680092373e-10 - 0.13494234691555555556e-12 * t) * t) * t) * t) * t) * t;
1822 }
1823 case 93: {
1824 double t = 2*y100 - 187;
1825 return 0.78191222288771379358e-1 + (-0.63914190297303976434e-2 + (0.29079759021299682675e-4 + (0.97885458059415717014e-8 + (-0.15635596116134296819e-8 + (0.17417110744051331974e-10 - 0.12591151763555555556e-12 * t) * t) * t) * t) * t) * t;
1826 }
1827 case 94: {
1828 double t = 2*y100 - 189;
1829 return 0.65524757106147402224e-1 + (-0.62750311956082444159e-2 + (0.29102328354323449795e-4 + (-0.20430838882727954582e-8 + (-0.13967781903855367270e-8 + (0.15958771833747057569e-10 - 0.11720175765333333333e-12 * t) * t) * t) * t) * t) * t;
1830 }
1831 case 95: {
1832 double t = 2*y100 - 191;
1833 return 0.53091065838453612773e-1 + (-0.61586898417077043662e-2 + (0.29057796072960100710e-4 + (-0.12597414620517987536e-7 + (-0.12440642607426861943e-8 + (0.14602787128447932137e-10 - 0.10885859114666666667e-12 * t) * t) * t) * t) * t) * t;
1834 }
1835 case 96: {
1836 double t = 2*y100 - 193;
1837 return 0.40889797115352738582e-1 + (-0.60426484889413678200e-2 + (0.28953496450191694606e-4 + (-0.21982952021823718400e-7 + (-0.11044169117553026211e-8 + (0.13344562332430552171e-10 - 0.10091231402844444444e-12 * t) * t) * t) * t) * t) * t;
1838 }
1839 case 97: case 98:
1840 case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
1841 // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9)
1842 double x2 = x*x;
1843 return x * (1.1283791670955125739
1844 - x2 * (0.75225277806367504925
1845 - x2 * (0.30090111122547001970
1846 - x2 * (0.085971746064420005629
1847 - x2 * 0.016931216931216931217))));
1848 }
1849 }
1850 /* Since 0 <= y100 < 101, this is only reached if x is NaN,
1851 in which case we should return NaN. */
1852 return NaN;
1853 }
1854
FADDEEVA(w_im)1855 double FADDEEVA(w_im)(double x)
1856 {
1857 if (x >= 0) {
1858 if (x > 45) { // continued-fraction expansion is faster
1859 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1860 if (x > 5e7) // 1-term expansion, important to avoid overflow
1861 return ispi / x;
1862 /* 5-term expansion (rely on compiler for CSE), simplified from:
1863 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1864 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1865 }
1866 return w_im_y100(100/(1+x), x);
1867 }
1868 else { // = -FADDEEVA(w_im)(-x)
1869 if (x < -45) { // continued-fraction expansion is faster
1870 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1871 if (x < -5e7) // 1-term expansion, important to avoid overflow
1872 return ispi / x;
1873 /* 5-term expansion (rely on compiler for CSE), simplified from:
1874 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1875 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1876 }
1877 return -w_im_y100(100/(1-x), -x);
1878 }
1879 }
1880
1881 /////////////////////////////////////////////////////////////////////////
1882
1883 // Compile with -DTEST_FADDEEVA to compile a little test program
1884 #if defined (TEST_FADDEEVA)
1885
1886 #if defined (__cplusplus)
1887 # include <cstdio>
1888 #else
1889 # include <stdio.h>
1890 #endif
1891
1892 // compute relative error |b-a|/|a|, handling case of NaN and Inf,
relerr(double a,double b)1893 static double relerr(double a, double b) {
1894 if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) {
1895 if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) ||
1896 (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) ||
1897 (isinf(a) && isinf(b) && a*b < 0))
1898 return Inf; // "infinite" error
1899 return 0; // matching infinity/nan results counted as zero error
1900 }
1901 if (a == 0)
1902 return b == 0 ? 0 : Inf;
1903 else
1904 return fabs((b-a) / a);
1905 }
1906
main(void)1907 int main(void) {
1908 double errmax_all = 0;
1909 {
1910 printf("############# w(z) tests #############\n");
1911 #define NTST 57 // define instead of const for C compatibility
1912 cmplx z[NTST] = {
1913 C(624.2,-0.26123),
1914 C(-0.4,3.),
1915 C(0.6,2.),
1916 C(-1.,1.),
1917 C(-1.,-9.),
1918 C(-1.,9.),
1919 C(-0.0000000234545,1.1234),
1920 C(-3.,5.1),
1921 C(-53,30.1),
1922 C(0.0,0.12345),
1923 C(11,1),
1924 C(-22,-2),
1925 C(9,-28),
1926 C(21,-33),
1927 C(1e5,1e5),
1928 C(1e14,1e14),
1929 C(-3001,-1000),
1930 C(1e160,-1e159),
1931 C(-6.01,0.01),
1932 C(-0.7,-0.7),
1933 C(2.611780000000000e+01, 4.540909610972489e+03),
1934 C(0.8e7,0.3e7),
1935 C(-20,-19.8081),
1936 C(1e-16,-1.1e-16),
1937 C(2.3e-8,1.3e-8),
1938 C(6.3,-1e-13),
1939 C(6.3,1e-20),
1940 C(1e-20,6.3),
1941 C(1e-20,16.3),
1942 C(9,1e-300),
1943 C(6.01,0.11),
1944 C(8.01,1.01e-10),
1945 C(28.01,1e-300),
1946 C(10.01,1e-200),
1947 C(10.01,-1e-200),
1948 C(10.01,0.99e-10),
1949 C(10.01,-0.99e-10),
1950 C(1e-20,7.01),
1951 C(-1,7.01),
1952 C(5.99,7.01),
1953 C(1,0),
1954 C(55,0),
1955 C(-0.1,0),
1956 C(1e-20,0),
1957 C(0,5e-14),
1958 C(0,51),
1959 C(Inf,0),
1960 C(-Inf,0),
1961 C(0,Inf),
1962 C(0,-Inf),
1963 C(Inf,Inf),
1964 C(Inf,-Inf),
1965 C(NaN,NaN),
1966 C(NaN,0),
1967 C(0,NaN),
1968 C(NaN,Inf),
1969 C(Inf,NaN)
1970 };
1971 cmplx w[NTST] = { /* w(z), computed with WolframAlpha
1972 ... note that WolframAlpha is problematic
1973 some of the above inputs, so I had to
1974 use the continued-fraction expansion
1975 in WolframAlpha in some cases, or switch
1976 to Maple */
1977 C(-3.78270245518980507452677445620103199303131110e-7,
1978 0.000903861276433172057331093754199933411710053155),
1979 C(0.1764906227004816847297495349730234591778719532788,
1980 -0.02146550539468457616788719893991501311573031095617),
1981 C(0.2410250715772692146133539023007113781272362309451,
1982 0.06087579663428089745895459735240964093522265589350),
1983 C(0.30474420525691259245713884106959496013413834051768,
1984 -0.20821893820283162728743734725471561394145872072738),
1985 C(7.317131068972378096865595229600561710140617977e34,
1986 8.321873499714402777186848353320412813066170427e34),
1987 C(0.0615698507236323685519612934241429530190806818395,
1988 -0.00676005783716575013073036218018565206070072304635),
1989 C(0.3960793007699874918961319170187598400134746631,
1990 -5.593152259116644920546186222529802777409274656e-9),
1991 C(0.08217199226739447943295069917990417630675021771804,
1992 -0.04701291087643609891018366143118110965272615832184),
1993 C(0.00457246000350281640952328010227885008541748668738,
1994 -0.00804900791411691821818731763401840373998654987934),
1995 C(0.8746342859608052666092782112565360755791467973338452,
1996 0.),
1997 C(0.00468190164965444174367477874864366058339647648741,
1998 0.0510735563901306197993676329845149741675029197050),
1999 C(-0.0023193175200187620902125853834909543869428763219,
2000 -0.025460054739731556004902057663500272721780776336),
2001 C(9.11463368405637174660562096516414499772662584e304,
2002 3.97101807145263333769664875189354358563218932e305),
2003 C(-4.4927207857715598976165541011143706155432296e281,
2004 -2.8019591213423077494444700357168707775769028e281),
2005 C(2.820947917809305132678577516325951485807107151e-6,
2006 2.820947917668257736791638444590253942253354058e-6),
2007 C(2.82094791773878143474039725787438662716372268e-15,
2008 2.82094791773878143474039725773333923127678361e-15),
2009 C(-0.0000563851289696244350147899376081488003110150498,
2010 -0.000169211755126812174631861529808288295454992688),
2011 C(-5.586035480670854326218608431294778077663867e-162,
2012 5.586035480670854326218608431294778077663867e-161),
2013 C(0.00016318325137140451888255634399123461580248456,
2014 -0.095232456573009287370728788146686162555021209999),
2015 C(0.69504753678406939989115375989939096800793577783885,
2016 -1.8916411171103639136680830887017670616339912024317),
2017 C(0.0001242418269653279656612334210746733213167234822,
2018 7.145975826320186888508563111992099992116786763e-7),
2019 C(2.318587329648353318615800865959225429377529825e-8,
2020 6.182899545728857485721417893323317843200933380e-8),
2021 C(-0.0133426877243506022053521927604277115767311800303,
2022 -0.0148087097143220769493341484176979826888871576145),
2023 C(1.00000000000000012412170838050638522857747934,
2024 1.12837916709551279389615890312156495593616433e-16),
2025 C(0.9999999853310704677583504063775310832036830015,
2026 2.595272024519678881897196435157270184030360773e-8),
2027 C(-1.4731421795638279504242963027196663601154624e-15,
2028 0.090727659684127365236479098488823462473074709),
2029 C(5.79246077884410284575834156425396800754409308e-18,
2030 0.0907276596841273652364790985059772809093822374),
2031 C(0.0884658993528521953466533278764830881245144368,
2032 1.37088352495749125283269718778582613192166760e-22),
2033 C(0.0345480845419190424370085249304184266813447878,
2034 2.11161102895179044968099038990446187626075258e-23),
2035 C(6.63967719958073440070225527042829242391918213e-36,
2036 0.0630820900592582863713653132559743161572639353),
2037 C(0.00179435233208702644891092397579091030658500743634,
2038 0.0951983814805270647939647438459699953990788064762),
2039 C(9.09760377102097999924241322094863528771095448e-13,
2040 0.0709979210725138550986782242355007611074966717),
2041 C(7.2049510279742166460047102593255688682910274423e-304,
2042 0.0201552956479526953866611812593266285000876784321),
2043 C(3.04543604652250734193622967873276113872279682e-44,
2044 0.0566481651760675042930042117726713294607499165),
2045 C(3.04543604652250734193622967873276113872279682e-44,
2046 0.0566481651760675042930042117726713294607499165),
2047 C(0.5659928732065273429286988428080855057102069081e-12,
2048 0.056648165176067504292998527162143030538756683302),
2049 C(-0.56599287320652734292869884280802459698927645e-12,
2050 0.0566481651760675042929985271621430305387566833029),
2051 C(0.0796884251721652215687859778119964009569455462,
2052 1.11474461817561675017794941973556302717225126e-22),
2053 C(0.07817195821247357458545539935996687005781943386550,
2054 -0.01093913670103576690766705513142246633056714279654),
2055 C(0.04670032980990449912809326141164730850466208439937,
2056 0.03944038961933534137558064191650437353429669886545),
2057 C(0.36787944117144232159552377016146086744581113103176,
2058 0.60715770584139372911503823580074492116122092866515),
2059 C(0,
2060 0.010259688805536830986089913987516716056946786526145),
2061 C(0.99004983374916805357390597718003655777207908125383,
2062 -0.11208866436449538036721343053869621153527769495574),
2063 C(0.99999999999999999999999999999999999999990000,
2064 1.12837916709551257389615890312154517168802603e-20),
2065 C(0.999999999999943581041645226871305192054749891144158,
2066 0),
2067 C(0.0110604154853277201542582159216317923453996211744250,
2068 0),
2069 C(0,0),
2070 C(0,0),
2071 C(0,0),
2072 C(Inf,0),
2073 C(0,0),
2074 C(NaN,NaN),
2075 C(NaN,NaN),
2076 C(NaN,NaN),
2077 C(NaN,0),
2078 C(NaN,NaN),
2079 C(NaN,NaN)
2080 };
2081 double errmax = 0;
2082 for (int i = 0; i < NTST; ++i) {
2083 cmplx fw = FADDEEVA(w)(z[i],0.);
2084 double re_err = relerr(creal(w[i]), creal(fw));
2085 double im_err = relerr(cimag(w[i]), cimag(fw));
2086 printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n",
2087 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]),
2088 re_err, im_err);
2089 if (re_err > errmax) errmax = re_err;
2090 if (im_err > errmax) errmax = im_err;
2091 }
2092 if (errmax > 1e-13) {
2093 printf("FAILURE -- relative error %g too large!\n", errmax);
2094 return 1;
2095 }
2096 printf("SUCCESS (max relative error = %g)\n", errmax);
2097 if (errmax > errmax_all) errmax_all = errmax;
2098 }
2099 {
2100 #undef NTST
2101 #define NTST 41 // define instead of const for C compatibility
2102 cmplx z[NTST] = {
2103 C(1,2),
2104 C(-1,2),
2105 C(1,-2),
2106 C(-1,-2),
2107 C(9,-28),
2108 C(21,-33),
2109 C(1e3,1e3),
2110 C(-3001,-1000),
2111 C(1e160,-1e159),
2112 C(5.1e-3, 1e-8),
2113 C(-4.9e-3, 4.95e-3),
2114 C(4.9e-3, 0.5),
2115 C(4.9e-4, -0.5e1),
2116 C(-4.9e-5, -0.5e2),
2117 C(5.1e-3, 0.5),
2118 C(5.1e-4, -0.5e1),
2119 C(-5.1e-5, -0.5e2),
2120 C(1e-6,2e-6),
2121 C(0,2e-6),
2122 C(0,2),
2123 C(0,20),
2124 C(0,200),
2125 C(Inf,0),
2126 C(-Inf,0),
2127 C(0,Inf),
2128 C(0,-Inf),
2129 C(Inf,Inf),
2130 C(Inf,-Inf),
2131 C(NaN,NaN),
2132 C(NaN,0),
2133 C(0,NaN),
2134 C(NaN,Inf),
2135 C(Inf,NaN),
2136 C(1e-3,NaN),
2137 C(7e-2,7e-2),
2138 C(7e-2,-7e-4),
2139 C(-9e-2,7e-4),
2140 C(-9e-2,9e-2),
2141 C(-7e-4,9e-2),
2142 C(7e-2,0.9e-2),
2143 C(7e-2,1.1e-2)
2144 };
2145 cmplx w[NTST] = { // erf(z[i]), evaluated with Maple
2146 C(-0.5366435657785650339917955593141927494421,
2147 -5.049143703447034669543036958614140565553),
2148 C(0.5366435657785650339917955593141927494421,
2149 -5.049143703447034669543036958614140565553),
2150 C(-0.5366435657785650339917955593141927494421,
2151 5.049143703447034669543036958614140565553),
2152 C(0.5366435657785650339917955593141927494421,
2153 5.049143703447034669543036958614140565553),
2154 C(0.3359473673830576996788000505817956637777e304,
2155 -0.1999896139679880888755589794455069208455e304),
2156 C(0.3584459971462946066523939204836760283645e278,
2157 0.3818954885257184373734213077678011282505e280),
2158 C(0.9996020422657148639102150147542224526887,
2159 0.00002801044116908227889681753993542916894856),
2160 C(-1, 0),
2161 C(1, 0),
2162 C(0.005754683859034800134412990541076554934877,
2163 0.1128349818335058741511924929801267822634e-7),
2164 C(-0.005529149142341821193633460286828381876955,
2165 0.005585388387864706679609092447916333443570),
2166 C(0.007099365669981359632319829148438283865814,
2167 0.6149347012854211635026981277569074001219),
2168 C(0.3981176338702323417718189922039863062440e8,
2169 -0.8298176341665249121085423917575122140650e10),
2170 C(-Inf,
2171 -Inf),
2172 C(0.007389128308257135427153919483147229573895,
2173 0.6149332524601658796226417164791221815139),
2174 C(0.4143671923267934479245651547534414976991e8,
2175 -0.8298168216818314211557046346850921446950e10),
2176 C(-Inf,
2177 -Inf),
2178 C(0.1128379167099649964175513742247082845155e-5,
2179 0.2256758334191777400570377193451519478895e-5),
2180 C(0,
2181 0.2256758334194034158904576117253481476197e-5),
2182 C(0,
2183 18.56480241457555259870429191324101719886),
2184 C(0,
2185 0.1474797539628786202447733153131835124599e173),
2186 C(0,
2187 Inf),
2188 C(1,0),
2189 C(-1,0),
2190 C(0,Inf),
2191 C(0,-Inf),
2192 C(NaN,NaN),
2193 C(NaN,NaN),
2194 C(NaN,NaN),
2195 C(NaN,0),
2196 C(0,NaN),
2197 C(NaN,NaN),
2198 C(NaN,NaN),
2199 C(NaN,NaN),
2200 C(0.07924380404615782687930591956705225541145,
2201 0.07872776218046681145537914954027729115247),
2202 C(0.07885775828512276968931773651224684454495,
2203 -0.0007860046704118224342390725280161272277506),
2204 C(-0.1012806432747198859687963080684978759881,
2205 0.0007834934747022035607566216654982820299469),
2206 C(-0.1020998418798097910247132140051062512527,
2207 0.1010030778892310851309082083238896270340),
2208 C(-0.0007962891763147907785684591823889484764272,
2209 0.1018289385936278171741809237435404896152),
2210 C(0.07886408666470478681566329888615410479530,
2211 0.01010604288780868961492224347707949372245),
2212 C(0.07886723099940260286824654364807981336591,
2213 0.01235199327873258197931147306290916629654)
2214 };
2215 #define TST(f,isc) \
2216 printf("############# " #f "(z) tests #############\n"); \
2217 double errmax = 0; \
2218 for (int i = 0; i < NTST; ++i) { \
2219 cmplx fw = FADDEEVA(f)(z[i],0.); \
2220 double re_err = relerr(creal(w[i]), creal(fw)); \
2221 double im_err = relerr(cimag(w[i]), cimag(fw)); \
2222 printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \
2223 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \
2224 re_err, im_err); \
2225 if (re_err > errmax) errmax = re_err; \
2226 if (im_err > errmax) errmax = im_err; \
2227 } \
2228 if (errmax > 1e-13) { \
2229 printf("FAILURE -- relative error %g too large!\n", errmax); \
2230 return 1; \
2231 } \
2232 printf("Checking " #f "(x) special case...\n"); \
2233 for (int i = 0; i < 10000; ++i) { \
2234 double x = pow(10., -300. + i * 600. / (10000 - 1)); \
2235 double re_err = relerr(FADDEEVA_RE(f)(x), \
2236 creal(FADDEEVA(f)(C(x,x*isc),0.))); \
2237 if (re_err > errmax) errmax = re_err; \
2238 re_err = relerr(FADDEEVA_RE(f)(-x), \
2239 creal(FADDEEVA(f)(C(-x,x*isc),0.))); \
2240 if (re_err > errmax) errmax = re_err; \
2241 } \
2242 { \
2243 double re_err = relerr(FADDEEVA_RE(f)(Inf), \
2244 creal(FADDEEVA(f)(C(Inf,0.),0.))); \
2245 if (re_err > errmax) errmax = re_err; \
2246 re_err = relerr(FADDEEVA_RE(f)(-Inf), \
2247 creal(FADDEEVA(f)(C(-Inf,0.),0.))); \
2248 if (re_err > errmax) errmax = re_err; \
2249 re_err = relerr(FADDEEVA_RE(f)(NaN), \
2250 creal(FADDEEVA(f)(C(NaN,0.),0.))); \
2251 if (re_err > errmax) errmax = re_err; \
2252 } \
2253 if (errmax > 1e-13) { \
2254 printf("FAILURE -- relative error %g too large!\n", errmax); \
2255 return 1; \
2256 } \
2257 printf("SUCCESS (max relative error = %g)\n", errmax); \
2258 if (errmax > errmax_all) errmax_all = errmax
2259
2260 TST(erf, 1e-20);
2261 }
2262 {
2263 // since erfi just calls through to erf, just one test should
2264 // be sufficient to make sure I didn't screw up the signs or something
2265 #undef NTST
2266 #define NTST 1 // define instead of const for C compatibility
2267 cmplx z[NTST] = { C(1.234,0.5678) };
2268 cmplx w[NTST] = { // erfi(z[i]), computed with Maple
2269 C(1.081032284405373149432716643834106923212,
2270 1.926775520840916645838949402886591180834)
2271 };
2272 TST(erfi, 0);
2273 }
2274 {
2275 // since erfcx just calls through to w, just one test should
2276 // be sufficient to make sure I didn't screw up the signs or something
2277 #undef NTST
2278 #define NTST 1 // define instead of const for C compatibility
2279 cmplx z[NTST] = { C(1.234,0.5678) };
2280 cmplx w[NTST] = { // erfcx(z[i]), computed with Maple
2281 C(0.3382187479799972294747793561190487832579,
2282 -0.1116077470811648467464927471872945833154)
2283 };
2284 TST(erfcx, 0);
2285 }
2286 {
2287 #undef NTST
2288 #define NTST 30 // define instead of const for C compatibility
2289 cmplx z[NTST] = {
2290 C(1,2),
2291 C(-1,2),
2292 C(1,-2),
2293 C(-1,-2),
2294 C(9,-28),
2295 C(21,-33),
2296 C(1e3,1e3),
2297 C(-3001,-1000),
2298 C(1e160,-1e159),
2299 C(5.1e-3, 1e-8),
2300 C(0,2e-6),
2301 C(0,2),
2302 C(0,20),
2303 C(0,200),
2304 C(2e-6,0),
2305 C(2,0),
2306 C(20,0),
2307 C(200,0),
2308 C(Inf,0),
2309 C(-Inf,0),
2310 C(0,Inf),
2311 C(0,-Inf),
2312 C(Inf,Inf),
2313 C(Inf,-Inf),
2314 C(NaN,NaN),
2315 C(NaN,0),
2316 C(0,NaN),
2317 C(NaN,Inf),
2318 C(Inf,NaN),
2319 C(88,0)
2320 };
2321 cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple
2322 C(1.536643565778565033991795559314192749442,
2323 5.049143703447034669543036958614140565553),
2324 C(0.4633564342214349660082044406858072505579,
2325 5.049143703447034669543036958614140565553),
2326 C(1.536643565778565033991795559314192749442,
2327 -5.049143703447034669543036958614140565553),
2328 C(0.4633564342214349660082044406858072505579,
2329 -5.049143703447034669543036958614140565553),
2330 C(-0.3359473673830576996788000505817956637777e304,
2331 0.1999896139679880888755589794455069208455e304),
2332 C(-0.3584459971462946066523939204836760283645e278,
2333 -0.3818954885257184373734213077678011282505e280),
2334 C(0.0003979577342851360897849852457775473112748,
2335 -0.00002801044116908227889681753993542916894856),
2336 C(2, 0),
2337 C(0, 0),
2338 C(0.9942453161409651998655870094589234450651,
2339 -0.1128349818335058741511924929801267822634e-7),
2340 C(1,
2341 -0.2256758334194034158904576117253481476197e-5),
2342 C(1,
2343 -18.56480241457555259870429191324101719886),
2344 C(1,
2345 -0.1474797539628786202447733153131835124599e173),
2346 C(1, -Inf),
2347 C(0.9999977432416658119838633199332831406314,
2348 0),
2349 C(0.004677734981047265837930743632747071389108,
2350 0),
2351 C(0.5395865611607900928934999167905345604088e-175,
2352 0),
2353 C(0, 0),
2354 C(0, 0),
2355 C(2, 0),
2356 C(1, -Inf),
2357 C(1, Inf),
2358 C(NaN, NaN),
2359 C(NaN, NaN),
2360 C(NaN, NaN),
2361 C(NaN, 0),
2362 C(1, NaN),
2363 C(NaN, NaN),
2364 C(NaN, NaN),
2365 C(0,0)
2366 };
2367 TST(erfc, 1e-20);
2368 }
2369 {
2370 #undef NTST
2371 #define NTST 48 // define instead of const for C compatibility
2372 cmplx z[NTST] = {
2373 C(2,1),
2374 C(-2,1),
2375 C(2,-1),
2376 C(-2,-1),
2377 C(-28,9),
2378 C(33,-21),
2379 C(1e3,1e3),
2380 C(-1000,-3001),
2381 C(1e-8, 5.1e-3),
2382 C(4.95e-3, -4.9e-3),
2383 C(5.1e-3, 5.1e-3),
2384 C(0.5, 4.9e-3),
2385 C(-0.5e1, 4.9e-4),
2386 C(-0.5e2, -4.9e-5),
2387 C(0.5e3, 4.9e-6),
2388 C(0.5, 5.1e-3),
2389 C(-0.5e1, 5.1e-4),
2390 C(-0.5e2, -5.1e-5),
2391 C(1e-6,2e-6),
2392 C(2e-6,0),
2393 C(2,0),
2394 C(20,0),
2395 C(200,0),
2396 C(0,4.9e-3),
2397 C(0,-5.1e-3),
2398 C(0,2e-6),
2399 C(0,-2),
2400 C(0,20),
2401 C(0,-200),
2402 C(Inf,0),
2403 C(-Inf,0),
2404 C(0,Inf),
2405 C(0,-Inf),
2406 C(Inf,Inf),
2407 C(Inf,-Inf),
2408 C(NaN,NaN),
2409 C(NaN,0),
2410 C(0,NaN),
2411 C(NaN,Inf),
2412 C(Inf,NaN),
2413 C(39, 6.4e-5),
2414 C(41, 6.09e-5),
2415 C(4.9e7, 5e-11),
2416 C(5.1e7, 4.8e-11),
2417 C(1e9, 2.4e-12),
2418 C(1e11, 2.4e-14),
2419 C(1e13, 2.4e-16),
2420 C(1e300, 2.4e-303)
2421 };
2422 cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple
2423 C(0.1635394094345355614904345232875688576839,
2424 -0.1531245755371229803585918112683241066853),
2425 C(-0.1635394094345355614904345232875688576839,
2426 -0.1531245755371229803585918112683241066853),
2427 C(0.1635394094345355614904345232875688576839,
2428 0.1531245755371229803585918112683241066853),
2429 C(-0.1635394094345355614904345232875688576839,
2430 0.1531245755371229803585918112683241066853),
2431 C(-0.01619082256681596362895875232699626384420,
2432 -0.005210224203359059109181555401330902819419),
2433 C(0.01078377080978103125464543240346760257008,
2434 0.006866888783433775382193630944275682670599),
2435 C(-0.5808616819196736225612296471081337245459,
2436 0.6688593905505562263387760667171706325749),
2437 C(Inf,
2438 -Inf),
2439 C(0.1000052020902036118082966385855563526705e-7,
2440 0.005100088434920073153418834680320146441685),
2441 C(0.004950156837581592745389973960217444687524,
2442 -0.004899838305155226382584756154100963570500),
2443 C(0.005100176864319675957314822982399286703798,
2444 0.005099823128319785355949825238269336481254),
2445 C(0.4244534840871830045021143490355372016428,
2446 0.002820278933186814021399602648373095266538),
2447 C(-0.1021340733271046543881236523269967674156,
2448 -0.00001045696456072005761498961861088944159916),
2449 C(-0.01000200120119206748855061636187197886859,
2450 0.9805885888237419500266621041508714123763e-8),
2451 C(0.001000002000012000023960527532953151819595,
2452 -0.9800058800588007290937355024646722133204e-11),
2453 C(0.4244549085628511778373438768121222815752,
2454 0.002935393851311701428647152230552122898291),
2455 C(-0.1021340732357117208743299813648493928105,
2456 -0.00001088377943049851799938998805451564893540),
2457 C(-0.01000200120119126652710792390331206563616,
2458 0.1020612612857282306892368985525393707486e-7),
2459 C(0.1000000000007333333333344266666666664457e-5,
2460 0.2000000000001333333333323199999999978819e-5),
2461 C(0.1999999999994666666666675199999999990248e-5,
2462 0),
2463 C(0.3013403889237919660346644392864226952119,
2464 0),
2465 C(0.02503136792640367194699495234782353186858,
2466 0),
2467 C(0.002500031251171948248596912483183760683918,
2468 0),
2469 C(0,0.004900078433419939164774792850907128053308),
2470 C(0,-0.005100088434920074173454208832365950009419),
2471 C(0,0.2000000000005333333333341866666666676419e-5),
2472 C(0,-48.16001211429122974789822893525016528191),
2473 C(0,0.4627407029504443513654142715903005954668e174),
2474 C(0,-Inf),
2475 C(0,0),
2476 C(-0,0),
2477 C(0, Inf),
2478 C(0, -Inf),
2479 C(NaN, NaN),
2480 C(NaN, NaN),
2481 C(NaN, NaN),
2482 C(NaN, 0),
2483 C(0, NaN),
2484 C(NaN, NaN),
2485 C(NaN, NaN),
2486 C(0.01282473148489433743567240624939698290584,
2487 -0.2105957276516618621447832572909153498104e-7),
2488 C(0.01219875253423634378984109995893708152885,
2489 -0.1813040560401824664088425926165834355953e-7),
2490 C(0.1020408163265306334945473399689037886997e-7,
2491 -0.1041232819658476285651490827866174985330e-25),
2492 C(0.9803921568627452865036825956835185367356e-8,
2493 -0.9227220299884665067601095648451913375754e-26),
2494 C(0.5000000000000000002500000000000000003750e-9,
2495 -0.1200000000000000001800000188712838420241e-29),
2496 C(5.00000000000000000000025000000000000000000003e-12,
2497 -1.20000000000000000000018000000000000000000004e-36),
2498 C(5.00000000000000000000000002500000000000000000e-14,
2499 -1.20000000000000000000000001800000000000000000e-42),
2500 C(5e-301, 0)
2501 };
2502 TST(Dawson, 1e-20);
2503 }
2504 printf("#####################################\n");
2505 printf("SUCCESS (max relative error = %g)\n", errmax_all);
2506 }
2507
2508 #endif
2509