File indexing completed on 2024-05-12 15:27:05
0001 /* Copyright (c) 2012 Massachusetts Institute of Technology 0002 * 0003 * Permission is hereby granted, free of charge, to any person obtaining 0004 * a copy of this software and associated documentation files (the 0005 * "Software"), to deal in the Software without restriction, including 0006 * without limitation the rights to use, copy, modify, merge, publish, 0007 * distribute, sublicense, and/or sell copies of the Software, and to 0008 * permit persons to whom the Software is furnished to do so, subject to 0009 * the following conditions: 0010 * 0011 * The above copyright notice and this permission notice shall be 0012 * included in all copies or substantial portions of the Software. 0013 * 0014 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 0015 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 0016 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 0017 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 0018 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 0019 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 0020 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 0021 */ 0022 0023 /* (Note that this file can be compiled with either C++, in which 0024 case it uses C++ std::complex<double>, or C, in which case it 0025 uses C99 double complex.) */ 0026 0027 /* Available at: http://ab-initio.mit.edu/Faddeeva 0028 0029 Computes various error functions (erf, erfc, erfi, erfcx), 0030 including the Dawson integral, in the complex plane, based 0031 on algorithms for the computation of the Faddeeva function 0032 w(z) = exp(-z^2) * erfc(-i*z). 0033 Given w(z), the error functions are mostly straightforward 0034 to compute, except for certain regions where we have to 0035 switch to Taylor expansions to avoid cancellation errors 0036 [e.g. near the origin for erf(z)]. 0037 0038 To compute the Faddeeva function, we use a combination of two 0039 algorithms: 0040 0041 For sufficiently large |z|, we use a continued-fraction expansion 0042 for w(z) similar to those described in: 0043 0044 Walter Gautschi, "Efficient computation of the complex error 0045 function," SIAM J. Numer. Anal. 7(1), pp. 187-198 (1970) 0046 0047 G. P. M. Poppe and C. M. J. Wijers, "More efficient computation 0048 of the complex error function," ACM Trans. Math. Soft. 16(1), 0049 pp. 38-46 (1990). 0050 0051 Unlike those papers, however, we switch to a completely different 0052 algorithm for smaller |z|: 0053 0054 Mofreh R. Zaghloul and Ahmed N. Ali, "Algorithm 916: Computing the 0055 Faddeyeva and Voigt Functions," ACM Trans. Math. Soft. 38(2), 15 0056 (2011). 0057 0058 (I initially used this algorithm for all z, but it turned out to be 0059 significantly slower than the continued-fraction expansion for 0060 larger |z|. On the other hand, it is competitive for smaller |z|, 0061 and is significantly more accurate than the Poppe & Wijers code 0062 in some regions, e.g. in the vicinity of z=1+1i.) 0063 0064 Note that this is an INDEPENDENT RE-IMPLEMENTATION of these algorithms, 0065 based on the description in the papers ONLY. In particular, I did 0066 not refer to the authors' Fortran or Matlab implementations, respectively, 0067 (which are under restrictive ACM copyright terms and therefore unusable 0068 in free/open-source software). 0069 0070 Steven G. Johnson, Massachusetts Institute of Technology 0071 https://math.mit.edu/~stevenj 0072 October 2012. 0073 0074 -- Note that Algorithm 916 assumes that the erfc(x) function, 0075 or rather the scaled function erfcx(x) = exp(x*x)*erfc(x), 0076 is supplied for REAL arguments x. I originally used an 0077 erfcx routine derived from DERFC in SLATEC, but I have 0078 since replaced it with a much faster routine written by 0079 me which uses a combination of continued-fraction expansions 0080 and a lookup table of Chebyshev polynomials. For speed, 0081 I implemented a similar algorithm for Im[w(x)] of real x, 0082 since this comes up frequently in the other error functions. 0083 0084 A small test program is included the end, which checks 0085 the w(z) etc. results against several known values. To compile 0086 the test function, compile with -DTEST_FADDEEVA (that is, 0087 #define TEST_FADDEEVA). 0088 0089 If HAVE_CONFIG_H is #defined (e.g. by compiling with -DHAVE_CONFIG_H), 0090 then we #include "config.h", which is assumed to be a GNU autoconf-style 0091 header defining HAVE_* macros to indicate the presence of features. In 0092 particular, if HAVE_ISNAN and HAVE_ISINF are #defined, we use those 0093 functions in math.h instead of defining our own, and if HAVE_ERF and/or 0094 HAVE_ERFC are defined we use those functions from <cmath> for erf and 0095 erfc of real arguments, respectively, instead of defining our own. 0096 0097 REVISION HISTORY: 0098 4 October 2012: Initial public release (SGJ) 0099 5 October 2012: Revised (SGJ) to fix spelling error, 0100 start summation for large x at round(x/a) (> 1) 0101 rather than ceil(x/a) as in the original 0102 paper, which should slightly improve performance 0103 (and, apparently, slightly improves accuracy) 0104 19 October 2012: Revised (SGJ) to fix bugs for large x, large -y, 0105 and 15<x<26. Performance improvements. Prototype 0106 now supplies default value for relerr. 0107 24 October 2012: Switch to continued-fraction expansion for 0108 sufficiently large z, for performance reasons. 0109 Also, avoid spurious overflow for |z| > 1e154. 0110 Set relerr argument to min(relerr,0.1). 0111 27 October 2012: Enhance accuracy in Re[w(z)] taken by itself, 0112 by switching to Alg. 916 in a region near 0113 the real-z axis where continued fractions 0114 have poor relative accuracy in Re[w(z)]. Thanks 0115 to M. Zaghloul for the tip. 0116 29 October 2012: Replace SLATEC-derived erfcx routine with 0117 completely rewritten code by me, using a very 0118 different algorithm which is much faster. 0119 30 October 2012: Implemented special-case code for real z 0120 (where real part is exp(-x^2) and imag part is 0121 Dawson integral), using algorithm similar to erfx. 0122 Export ImFaddeeva_w function to make Dawson's 0123 integral directly accessible. 0124 3 November 2012: Provide implementations of erf, erfc, erfcx, 0125 and Dawson functions in Faddeeva:: namespace, 0126 in addition to Faddeeva::w. Provide header 0127 file Faddeeva.hh. 0128 4 November 2012: Slightly faster erf for real arguments. 0129 Updated MATLAB and Octave plugins. 0130 27 November 2012: Support compilation with either C++ or 0131 plain C (using C99 complex numbers). 0132 For real x, use standard-library erf(x) 0133 and erfc(x) if available (for C99 or C++11). 0134 #include "config.h" if HAVE_CONFIG_H is #defined. 0135 15 December 2012: Portability fixes (copysign, Inf/NaN creation), 0136 use CMPLX/__builtin_complex if available in C, 0137 slight accuracy improvements to erf and dawson 0138 functions near the origin. Use gnulib functions 0139 if GNULIB_NAMESPACE is defined. 0140 18 December 2012: Slight tweaks (remove recomputation of x*x in Dawson) 0141 12 May 2015: Bugfix for systems lacking copysign function. 0142 */ 0143 0144 ///////////////////////////////////////////////////////////////////////// 0145 /* If this file is compiled as a part of a larger project, 0146 support using an autoconf-style config.h header file 0147 (with various "HAVE_*" #defines to indicate features) 0148 if HAVE_CONFIG_H is #defined (in GNU autotools style). */ 0149 0150 #ifdef HAVE_CONFIG_H 0151 # include "config.h" 0152 #endif 0153 0154 ///////////////////////////////////////////////////////////////////////// 0155 // macros to allow us to use either C++ or C (with C99 features) 0156 0157 #ifdef __cplusplus 0158 0159 # include "Faddeeva.hh" 0160 0161 # include <cfloat> 0162 # include <cmath> 0163 # include <limits> 0164 using namespace std; 0165 0166 // use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS) 0167 # define Inf numeric_limits<double>::infinity() 0168 # define NaN numeric_limits<double>::quiet_NaN() 0169 0170 typedef complex<double> cmplx; 0171 0172 // Use C-like complex syntax, since the C syntax is more restrictive 0173 # define cexp(z) exp(z) 0174 # define creal(z) real(z) 0175 # define cimag(z) imag(z) 0176 # define cpolar(r,t) polar(r,t) 0177 0178 # define C(a,b) cmplx(a,b) 0179 0180 # define FADDEEVA(name) Faddeeva::name 0181 # define FADDEEVA_RE(name) Faddeeva::name 0182 0183 // isnan/isinf were introduced in C++11 0184 # if (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF)) 0185 static inline bool my_isnan(double x) { return x != x; } 0186 # define isnan my_isnan 0187 static inline bool my_isinf(double x) { return 1/x == 0.; } 0188 # define isinf my_isinf 0189 # elif (__cplusplus >= 201103L) 0190 // g++ gets confused between the C and C++ isnan/isinf functions 0191 # define isnan std::isnan 0192 # define isinf std::isinf 0193 # endif 0194 0195 // copysign was introduced in C++11 (and is also in POSIX and C99) 0196 # if defined(_WIN32) || defined(__WIN32__) 0197 # define copysign _copysign // of course MS had to be different 0198 # elif defined(GNULIB_NAMESPACE) // we are using using gnulib <cmath> 0199 # define copysign GNULIB_NAMESPACE::copysign 0200 # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX) 0201 static inline double my_copysign(double x, double y) { return x<0 != y<0 ? -x : x; } 0202 # define copysign my_copysign 0203 # endif 0204 0205 // If we are using the gnulib <cmath> (e.g. in the GNU Octave sources), 0206 // gnulib generates a link warning if we use ::floor instead of gnulib::floor. 0207 // This warning is completely innocuous because the only difference between 0208 // gnulib::floor and the system ::floor (and only on ancient OSF systems) 0209 // has to do with floor(-0), which doesn't occur in the usage below, but 0210 // the Octave developers prefer that we silence the warning. 0211 # ifdef GNULIB_NAMESPACE 0212 # define floor GNULIB_NAMESPACE::floor 0213 # endif 0214 0215 #else // !__cplusplus, i.e. pure C (requires C99 features) 0216 0217 # include "Faddeeva.h" 0218 0219 # ifndef _GNU_SOURCE 0220 # define _GNU_SOURCE // enable GNU libc NAN extension if possible 0221 # endif 0222 0223 # include <float.h> 0224 # include <math.h> 0225 0226 typedef double complex cmplx; 0227 0228 # define FADDEEVA(name) Faddeeva_ ## name 0229 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re 0230 0231 /* Constructing complex numbers like 0+i*NaN is problematic in C99 0232 without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if 0233 I is a complex (rather than imaginary) constant. For some reason, 0234 however, it works fine in (pre-4.7) gcc if I define Inf and NaN as 0235 1/0 and 0/0 (and only if I compile with optimization -O1 or more), 0236 but not if I use the INFINITY or NAN macros. */ 0237 0238 /* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro 0239 may not be defined unless we are using a recent (2012) version of 0240 glibc and compile with -std=c11... note that icc lies about being 0241 gcc and probably doesn't have this builtin(?), so exclude icc explicitly */ 0242 # if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER)) 0243 # define CMPLX(a,b) __builtin_complex((double) (a), (double) (b)) 0244 # endif 0245 0246 # ifdef CMPLX // C11 0247 # define C(a,b) CMPLX(a,b) 0248 # define Inf INFINITY // C99 infinity 0249 # ifdef NAN // GNU libc extension 0250 # define NaN NAN 0251 # else 0252 # define NaN (0./0.) // NaN 0253 # endif 0254 # else 0255 # define C(a,b) ((a) + I*(b)) 0256 # define Inf (1./0.) 0257 # define NaN (0./0.) 0258 # endif 0259 0260 static inline cmplx cpolar(double r, double t) 0261 { 0262 if (r == 0.0 && !isnan(t)) 0263 return 0.0; 0264 else 0265 return C(r * cos(t), r * sin(t)); 0266 } 0267 0268 #endif // !__cplusplus, i.e. pure C (requires C99 features) 0269 0270 ///////////////////////////////////////////////////////////////////////// 0271 // Auxiliary routines to compute other special functions based on w(z) 0272 0273 // compute erfcx(z) = exp(z^2) erfz(z) 0274 cmplx FADDEEVA(erfcx)(cmplx z, double relerr) 0275 { 0276 return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr); 0277 } 0278 0279 // compute the error function erf(x) 0280 double FADDEEVA_RE(erf)(double x) 0281 { 0282 #if !defined(__cplusplus) 0283 return erf(x); // C99 supplies erf in math.h 0284 #elif (defined(__cplusplus) && __cplusplus >= 201103L) || defined(HAVE_ERF) 0285 return ::erf(x); // C++11 supplies std::erf in cmath 0286 #else 0287 double mx2 = -x*x; 0288 if (mx2 < -750) // underflow 0289 return (x >= 0 ? 1.0 : -1.0); 0290 0291 if (x >= 0) { 0292 if (x < 8e-2) goto taylor; 0293 return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x); 0294 } 0295 else { // x < 0 0296 if (x > -8e-2) goto taylor; 0297 return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0; 0298 } 0299 0300 // Use Taylor series for small |x|, to avoid cancellation inaccuracy 0301 // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...) 0302 taylor: 0303 return x * (1.1283791670955125739 0304 + mx2 * (0.37612638903183752464 0305 + mx2 * (0.11283791670955125739 0306 + mx2 * (0.026866170645131251760 0307 + mx2 * 0.0052239776254421878422)))); 0308 #endif 0309 } 0310 0311 // compute the error function erf(z) 0312 cmplx FADDEEVA(erf)(cmplx z, double relerr) 0313 { 0314 double x = creal(z), y = cimag(z); 0315 0316 if (y == 0) 0317 return C(FADDEEVA_RE(erf)(x), 0318 y); // preserve sign of 0 0319 if (x == 0) // handle separately for speed & handling of y = Inf or NaN 0320 return C(x, // preserve sign of 0 0321 /* handle y -> Inf limit manually, since 0322 exp(y^2) -> Inf but Im[w(y)] -> 0, so 0323 IEEE will give us a NaN when it should be Inf */ 0324 y*y > 720 ? (y > 0 ? Inf : -Inf) 0325 : exp(y*y) * FADDEEVA(w_im)(y)); 0326 0327 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 0328 double mIm_z2 = -2*x*y; // Im(-z^2) 0329 if (mRe_z2 < -750) // underflow 0330 return (x >= 0 ? 1.0 : -1.0); 0331 0332 /* Handle positive and negative x via different formulas, 0333 using the mirror symmetries of w, to avoid overflow/underflow 0334 problems from multiplying exponentially large and small quantities. */ 0335 if (x >= 0) { 0336 if (x < 8e-2) { 0337 if (fabs(y) < 1e-2) 0338 goto taylor; 0339 else if (fabs(mIm_z2) < 5e-3 && x < 5e-3) 0340 goto taylor_erfi; 0341 } 0342 /* don't use complex exp function, since that will produce spurious NaN 0343 values when multiplying w in an overflow situation. */ 0344 return 1.0 - exp(mRe_z2) * 0345 (C(cos(mIm_z2), sin(mIm_z2)) 0346 * FADDEEVA(w)(C(-y,x), relerr)); 0347 } 0348 else { // x < 0 0349 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call 0350 if (fabs(y) < 1e-2) 0351 goto taylor; 0352 else if (fabs(mIm_z2) < 5e-3 && x > -5e-3) 0353 goto taylor_erfi; 0354 } 0355 else if (isnan(x)) 0356 return C(NaN, y == 0 ? 0 : NaN); 0357 /* don't use complex exp function, since that will produce spurious NaN 0358 values when multiplying w in an overflow situation. */ 0359 return exp(mRe_z2) * 0360 (C(cos(mIm_z2), sin(mIm_z2)) 0361 * FADDEEVA(w)(C(y,-x), relerr)) - 1.0; 0362 } 0363 0364 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 0365 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...) 0366 taylor: 0367 { 0368 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 0369 return z * (1.1283791670955125739 0370 + mz2 * (0.37612638903183752464 0371 + mz2 * (0.11283791670955125739 0372 + mz2 * (0.026866170645131251760 0373 + mz2 * 0.0052239776254421878422)))); 0374 } 0375 0376 /* for small |x| and small |xy|, 0377 use Taylor series to avoid cancellation inaccuracy: 0378 erf(x+iy) = erf(iy) 0379 + 2*exp(y^2)/sqrt(pi) * 0380 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 0381 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] 0382 where: 0383 erf(iy) = exp(y^2) * Im[w(y)] 0384 */ 0385 taylor_erfi: 0386 { 0387 double x2 = x*x, y2 = y*y; 0388 double expy2 = exp(y2); 0389 return C 0390 (expy2 * x * (1.1283791670955125739 0391 - x2 * (0.37612638903183752464 0392 + 0.75225277806367504925*y2) 0393 + x2*x2 * (0.11283791670955125739 0394 + y2 * (0.45135166683820502956 0395 + 0.15045055561273500986*y2))), 0396 expy2 * (FADDEEVA(w_im)(y) 0397 - x2*y * (1.1283791670955125739 0398 - x2 * (0.56418958354775628695 0399 + 0.37612638903183752464*y2)))); 0400 } 0401 } 0402 0403 // erfi(z) = -i erf(iz) 0404 cmplx FADDEEVA(erfi)(cmplx z, double relerr) 0405 { 0406 cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr); 0407 return C(cimag(e), -creal(e)); 0408 } 0409 0410 // erfi(x) = -i erf(ix) 0411 double FADDEEVA_RE(erfi)(double x) 0412 { 0413 return x*x > 720 ? (x > 0 ? Inf : -Inf) 0414 : exp(x*x) * FADDEEVA(w_im)(x); 0415 } 0416 0417 // erfc(x) = 1 - erf(x) 0418 double FADDEEVA_RE(erfc)(double x) 0419 { 0420 #if !defined(__cplusplus) 0421 return erfc(x); // C99 supplies erfc in math.h 0422 #elif (defined(__cplusplus) && __cplusplus >= 201103L) || defined(HAVE_ERFC) 0423 return ::erfc(x); // C++11 supplies std::erfc in cmath 0424 #else 0425 if (x*x > 750) // underflow 0426 return (x >= 0 ? 0.0 : 2.0); 0427 return x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 0428 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x); 0429 #endif 0430 } 0431 0432 // erfc(z) = 1 - erf(z) 0433 cmplx FADDEEVA(erfc)(cmplx z, double relerr) 0434 { 0435 double x = creal(z), y = cimag(z); 0436 0437 if (x == 0.) 0438 return C(1, 0439 /* handle y -> Inf limit manually, since 0440 exp(y^2) -> Inf but Im[w(y)] -> 0, so 0441 IEEE will give us a NaN when it should be Inf */ 0442 y*y > 720 ? (y > 0 ? -Inf : Inf) 0443 : -exp(y*y) * FADDEEVA(w_im)(y)); 0444 if (y == 0.) { 0445 if (x*x > 750) // underflow 0446 return C(x >= 0 ? 0.0 : 2.0, 0447 -y); // preserve sign of 0 0448 return C(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 0449 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x), 0450 -y); // preserve sign of zero 0451 } 0452 0453 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 0454 double mIm_z2 = -2*x*y; // Im(-z^2) 0455 if (mRe_z2 < -750) // underflow 0456 return (x >= 0 ? 0.0 : 2.0); 0457 0458 if (x >= 0) 0459 return cexp(C(mRe_z2, mIm_z2)) 0460 * FADDEEVA(w)(C(-y,x), relerr); 0461 else 0462 return 2.0 - cexp(C(mRe_z2, mIm_z2)) 0463 * FADDEEVA(w)(C(y,-x), relerr); 0464 } 0465 0466 // compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x) 0467 double FADDEEVA_RE(Dawson)(double x) 0468 { 0469 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 0470 return spi2 * FADDEEVA(w_im)(x); 0471 } 0472 0473 // compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z) 0474 cmplx FADDEEVA(Dawson)(cmplx z, double relerr) 0475 { 0476 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 0477 double x = creal(z), y = cimag(z); 0478 0479 // handle axes separately for speed & proper handling of x or y = Inf or NaN 0480 if (y == 0) 0481 return C(spi2 * FADDEEVA(w_im)(x), 0482 -y); // preserve sign of 0 0483 if (x == 0) { 0484 double y2 = y*y; 0485 if (y2 < 2.5e-5) { // Taylor expansion 0486 return C(x, // preserve sign of 0 0487 y * (1. 0488 + y2 * (0.6666666666666666666666666666666666666667 0489 + y2 * 0.26666666666666666666666666666666666667))); 0490 } 0491 return C(x, // preserve sign of 0 0492 spi2 * (y >= 0 0493 ? exp(y2) - FADDEEVA_RE(erfcx)(y) 0494 : FADDEEVA_RE(erfcx)(-y) - exp(y2))); 0495 } 0496 0497 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 0498 double mIm_z2 = -2*x*y; // Im(-z^2) 0499 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 0500 0501 /* Handle positive and negative x via different formulas, 0502 using the mirror symmetries of w, to avoid overflow/underflow 0503 problems from multiplying exponentially large and small quantities. */ 0504 if (y >= 0) { 0505 if (y < 5e-3) { 0506 if (fabs(x) < 5e-3) 0507 goto taylor; 0508 else if (fabs(mIm_z2) < 5e-3) 0509 goto taylor_realaxis; 0510 } 0511 cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr); 0512 return spi2 * C(-cimag(res), creal(res)); 0513 } 0514 else { // y < 0 0515 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call 0516 if (fabs(x) < 5e-3) 0517 goto taylor; 0518 else if (fabs(mIm_z2) < 5e-3) 0519 goto taylor_realaxis; 0520 } 0521 else if (isnan(y)) 0522 return C(x == 0 ? 0 : NaN, NaN); 0523 cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2); 0524 return spi2 * C(-cimag(res), creal(res)); 0525 } 0526 0527 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 0528 // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... 0529 taylor: 0530 return z * (1. 0531 + mz2 * (0.6666666666666666666666666666666666666667 0532 + mz2 * 0.2666666666666666666666666666666666666667)); 0533 0534 /* for small |y| and small |xy|, 0535 use Taylor series to avoid cancellation inaccuracy: 0536 dawson(x + iy) 0537 = D + y^2 (D + x - 2Dx^2) 0538 + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3) 0539 + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3) 0540 + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ... 0541 where D = dawson(x) 0542 0543 However, for large |x|, 2Dx -> 1 which gives cancellation problems in 0544 this series (many of the leading terms cancel). So, for large |x|, 0545 we need to substitute a continued-fraction expansion for D. 0546 0547 dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...)))))) 0548 0549 The 6 terms shown here seems to be the minimum needed to be 0550 accurate as soon as the simpler Taylor expansion above starts 0551 breaking down. Using this 6-term expansion, factoring out the 0552 denominator, and simplifying with Maple, we obtain: 0553 0554 Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x 0555 = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4 0556 Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y 0557 = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4 0558 0559 Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction 0560 expansion for the real part, and a 2-term expansion for the imaginary 0561 part. (This avoids overflow problems for huge |x|.) This yields: 0562 0563 Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x) 0564 Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1) 0565 0566 */ 0567 taylor_realaxis: 0568 { 0569 double x2 = x*x; 0570 if (x2 > 1600) { // |x| > 40 0571 double y2 = y*y; 0572 if (x2 > 25e14) {// |x| > 5e7 0573 double xy2 = (x*y)*(x*y); 0574 return C((0.5 + y2 * (0.5 + 0.25*y2 0575 - 0.16666666666666666667*xy2)) / x, 0576 y * (-1 + y2 * (-0.66666666666666666667 0577 + 0.13333333333333333333*xy2 0578 - 0.26666666666666666667*y2)) 0579 / (2*x2 - 1)); 0580 } 0581 return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) * 0582 C(x * (33 + x2 * (-28 + 4*x2) 0583 + y2 * (18 - 4*x2 + 4*y2)), 0584 y * (-15 + x2 * (24 - 4*x2) 0585 + y2 * (4*x2 - 10 - 4*y2))); 0586 } 0587 else { 0588 double D = spi2 * FADDEEVA(w_im)(x); 0589 double y2 = y*y; 0590 return C 0591 (D + y2 * (D + x - 2*D*x2) 0592 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2)) 0593 + x * (0.83333333333333333333 0594 - 0.33333333333333333333 * x2)), 0595 y * (1 - 2*D*x 0596 + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2)) 0597 + y2*y2 * (0.26666666666666666667 - 0598 x2 * (0.6 - 0.13333333333333333333 * x2) 0599 - D*x * (1 - x2 * (1.3333333333333333333 0600 - 0.26666666666666666667 * x2))))); 0601 } 0602 } 0603 } 0604 0605 ///////////////////////////////////////////////////////////////////////// 0606 0607 // return sinc(x) = sin(x)/x, given both x and sin(x) 0608 // [since we only use this in cases where sin(x) has already been computed] 0609 static inline double sinc(double x, double sinx) { 0610 return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x; 0611 } 0612 0613 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2 0614 static inline double sinh_taylor(double x) { 0615 return x * (1 + (x*x) * (0.1666666666666666666667 0616 + 0.00833333333333333333333 * (x*x))); 0617 } 0618 0619 static inline double sqr(double x) { return x*x; } 0620 0621 // precomputed table of expa2n2[n-1] = exp(-a2*n*n) 0622 // for double-precision a2 = 0.26865... in FADDEEVA(w), below. 0623 const int expa2n2_size=52; 0624 static const double expa2n2[] = { 0625 7.64405281671221563e-01, 0626 3.41424527166548425e-01, 0627 8.91072646929412548e-02, 0628 1.35887299055460086e-02, 0629 1.21085455253437481e-03, 0630 6.30452613933449404e-05, 0631 1.91805156577114683e-06, 0632 3.40969447714832381e-08, 0633 3.54175089099469393e-10, 0634 2.14965079583260682e-12, 0635 7.62368911833724354e-15, 0636 1.57982797110681093e-17, 0637 1.91294189103582677e-20, 0638 1.35344656764205340e-23, 0639 5.59535712428588720e-27, 0640 1.35164257972401769e-30, 0641 1.90784582843501167e-34, 0642 1.57351920291442930e-38, 0643 7.58312432328032845e-43, 0644 2.13536275438697082e-47, 0645 3.51352063787195769e-52, 0646 3.37800830266396920e-57, 0647 1.89769439468301000e-62, 0648 6.22929926072668851e-68, 0649 1.19481172006938722e-73, 0650 1.33908181133005953e-79, 0651 8.76924303483223939e-86, 0652 3.35555576166254986e-92, 0653 7.50264110688173024e-99, 0654 9.80192200745410268e-106, 0655 7.48265412822268959e-113, 0656 3.33770122566809425e-120, 0657 8.69934598159861140e-128, 0658 1.32486951484088852e-135, 0659 1.17898144201315253e-143, 0660 6.13039120236180012e-152, 0661 1.86258785950822098e-160, 0662 3.30668408201432783e-169, 0663 3.43017280887946235e-178, 0664 2.07915397775808219e-187, 0665 7.36384545323984966e-197, 0666 1.52394760394085741e-206, 0667 1.84281935046532100e-216, 0668 1.30209553802992923e-226, 0669 5.37588903521080531e-237, 0670 1.29689584599763145e-247, 0671 1.82813078022866562e-258, 0672 1.50576355348684241e-269, 0673 7.24692320799294194e-281, 0674 2.03797051314726829e-292, 0675 3.34880215927873807e-304, 0676 0.0 // underflow (also prevents reads past array end, below) 0677 }; 0678 0679 ///////////////////////////////////////////////////////////////////////// 0680 0681 cmplx FADDEEVA(w)(cmplx z, double relerr) 0682 { 0683 if (creal(z) == 0.0) 0684 return C(FADDEEVA_RE(erfcx)(cimag(z)), 0685 creal(z)); // give correct sign of 0 in cimag(w) 0686 else if (cimag(z) == 0) 0687 return C(exp(-sqr(creal(z))), 0688 FADDEEVA(w_im)(creal(z))); 0689 0690 double a, a2, c; 0691 if (relerr <= DBL_EPSILON) { 0692 relerr = DBL_EPSILON; 0693 a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5)) 0694 c = 0.329973702884629072537; // (2/pi) * a; 0695 a2 = 0.268657157075235951582; // a^2 0696 } 0697 else { 0698 const double pi = 3.14159265358979323846264338327950288419716939937510582; 0699 if (relerr > 0.1) relerr = 0.1; // not sensible to compute < 1 digit 0700 a = pi / sqrt(-log(relerr*0.5)); 0701 c = (2/pi)*a; 0702 a2 = a*a; 0703 } 0704 const double x = fabs(creal(z)); 0705 const double y = cimag(z), ya = fabs(y); 0706 0707 cmplx ret = 0.; // return value 0708 0709 double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0; 0710 0711 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z| 0712 0713 #if USE_CONTINUED_FRACTION 0714 if (ya > 7 || (x > 6 // continued fraction is faster 0715 /* As pointed out by M. Zaghloul, the continued 0716 fraction seems to give a large relative error in 0717 Re w(z) for |x| ~ 6 and small |y|, so use 0718 algorithm 816 in this region: */ 0719 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) { 0720 0721 /* Poppe & Wijers suggest using a number of terms 0722 nu = 3 + 1442 / (26*rho + 77) 0723 where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4. 0724 (They only use this expansion for rho >= 1, but rho a little less 0725 than 1 seems okay too.) 0726 Instead, I did my own fit to a slightly different function 0727 that avoids the hypotenuse calculation, using NLopt to minimize 0728 the sum of the squares of the errors in nu with the constraint 0729 that the estimated nu be >= minimum nu to attain machine precision. 0730 I also separate the regions where nu == 2 and nu == 1. */ 0731 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 0732 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 0733 if (x + ya > 4000) { // nu <= 2 0734 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z 0735 // scale to avoid overflow 0736 if (x > ya) { 0737 double yax = ya / xs; 0738 double denom = ispi / (xs + yax*ya); 0739 ret = C(denom*yax, denom); 0740 } 0741 else if (isinf(ya)) 0742 return ((isnan(x) || y < 0) 0743 ? C(NaN,NaN) : C(0,0)); 0744 else { 0745 double xya = xs / ya; 0746 double denom = ispi / (xya*xs + ya); 0747 ret = C(denom, denom*xya); 0748 } 0749 } 0750 else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5) 0751 double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya; 0752 double denom = ispi / (dr*dr + di*di); 0753 ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di)); 0754 } 0755 } 0756 else { // compute nu(z) estimate and do general continued fraction 0757 const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit 0758 double nu = floor(c0 + c1 / (c2*x + c3*ya + c4)); 0759 double wr = xs, wi = ya; 0760 for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) { 0761 // w <- z - nu/w: 0762 double denom = nu / (wr*wr + wi*wi); 0763 wr = xs - wr * denom; 0764 wi = ya + wi * denom; 0765 } 0766 { // w(z) = i/sqrt(pi) / w: 0767 double denom = ispi / (wr*wr + wi*wi); 0768 ret = C(denom*wi, denom*wr); 0769 } 0770 } 0771 if (y < 0) { 0772 // use w(z) = 2.0*exp(-z*z) - w(-z), 0773 // but be careful of overflow in exp(-z*z) 0774 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) 0775 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret; 0776 } 0777 else 0778 return ret; 0779 } 0780 #else // !USE_CONTINUED_FRACTION 0781 if (x + ya > 1e7) { // w(z) = i/sqrt(pi) / z, to machine precision 0782 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 0783 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 0784 // scale to avoid overflow 0785 if (x > ya) { 0786 double yax = ya / xs; 0787 double denom = ispi / (xs + yax*ya); 0788 ret = C(denom*yax, denom); 0789 } 0790 else { 0791 double xya = xs / ya; 0792 double denom = ispi / (xya*xs + ya); 0793 ret = C(denom, denom*xya); 0794 } 0795 if (y < 0) { 0796 // use w(z) = 2.0*exp(-z*z) - w(-z), 0797 // but be careful of overflow in exp(-z*z) 0798 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) 0799 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret; 0800 } 0801 else 0802 return ret; 0803 } 0804 #endif // !USE_CONTINUED_FRACTION 0805 0806 /* Note: The test that seems to be suggested in the paper is x < 0807 sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2) 0808 underflows to zero and sum1,sum2,sum4 are zero. However, long 0809 before this occurs, the sum1,sum2,sum4 contributions are 0810 negligible in double precision; I find that this happens for x > 0811 about 6, for all y. On the other hand, I find that the case 0812 where we compute all of the sums is faster (at least with the 0813 precomputed expa2n2 table) until about x=10. Furthermore, if we 0814 try to compute all of the sums for x > 20, I find that we 0815 sometimes run into numerical problems because underflow/overflow 0816 problems start to appear in the various coefficients of the sums, 0817 below. Therefore, we use x < 10 here. */ 0818 else if (x < 10) { 0819 double prod2ax = 1, prodm2ax = 1; 0820 double expx2; 0821 0822 if (isnan(y)) 0823 return C(y,y); 0824 0825 /* Somewhat ugly copy-and-paste duplication here, but I see significant 0826 speedups from using the special-case code with the precomputed 0827 exponential, and the x < 5e-4 special case is needed for accuracy. */ 0828 0829 if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table 0830 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 0831 const double x2 = x*x; 0832 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor 0833 // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision 0834 const double ax2 = 1.036642960860171859744*x; // 2*a*x 0835 const double exp2ax = 0836 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2)); 0837 const double expm2ax = 0838 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2)); 0839 for (int n = 1; n <= expa2n2_size; ++n) { 0840 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); 0841 prod2ax *= exp2ax; 0842 prodm2ax *= expm2ax; 0843 sum1 += coef; 0844 sum2 += coef * prodm2ax; 0845 sum3 += coef * prod2ax; 0846 0847 // really = sum5 - sum4 0848 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x); 0849 0850 // test convergence via sum3 0851 if (coef * prod2ax < relerr * sum3) break; 0852 } 0853 } 0854 else { // x > 5e-4, compute sum4 and sum5 separately 0855 expx2 = exp(-x*x); 0856 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax; 0857 for (int n = 1; n <= expa2n2_size; ++n) { 0858 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); 0859 prod2ax *= exp2ax; 0860 prodm2ax *= expm2ax; 0861 sum1 += coef; 0862 sum2 += coef * prodm2ax; 0863 sum4 += (coef * prodm2ax) * (a*n); 0864 sum3 += coef * prod2ax; 0865 sum5 += (coef * prod2ax) * (a*n); 0866 // test convergence via sum5, since this sum has the slowest decay 0867 if ((coef * prod2ax) * (a*n) < relerr * sum5) break; 0868 } 0869 } 0870 } 0871 else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly 0872 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax; 0873 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 0874 const double x2 = x*x; 0875 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor 0876 for (int n = 1; 1; ++n) { 0877 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y); 0878 prod2ax *= exp2ax; 0879 prodm2ax *= expm2ax; 0880 sum1 += coef; 0881 sum2 += coef * prodm2ax; 0882 sum3 += coef * prod2ax; 0883 0884 // really = sum5 - sum4 0885 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x); 0886 0887 // test convergence via sum3 0888 if (coef * prod2ax < relerr * sum3) break; 0889 } 0890 } 0891 else { // x > 5e-4, compute sum4 and sum5 separately 0892 expx2 = exp(-x*x); 0893 for (int n = 1; 1; ++n) { 0894 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y); 0895 prod2ax *= exp2ax; 0896 prodm2ax *= expm2ax; 0897 sum1 += coef; 0898 sum2 += coef * prodm2ax; 0899 sum4 += (coef * prodm2ax) * (a*n); 0900 sum3 += coef * prod2ax; 0901 sum5 += (coef * prod2ax) * (a*n); 0902 // test convergence via sum5, since this sum has the slowest decay 0903 if ((coef * prod2ax) * (a*n) < relerr * sum5) break; 0904 } 0905 } 0906 } 0907 const double expx2erfcxy = // avoid spurious overflow for large negative y 0908 y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision 0909 ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x); 0910 if (y > 5) { // imaginary terms cancel 0911 const double sinxy = sin(x*y); 0912 ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y) 0913 + (c*x*expx2) * sinxy * sinc(x*y, sinxy); 0914 } 0915 else { 0916 double xs = creal(z); 0917 const double sinxy = sin(xs*y); 0918 const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y); 0919 const double coef1 = expx2erfcxy - c*y*sum1; 0920 const double coef2 = c*xs*expx2; 0921 ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy), 0922 coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy); 0923 } 0924 } 0925 else { // x large: only sum3 & sum5 contribute (see above note) 0926 if (isnan(x)) 0927 return C(x,x); 0928 if (isnan(y)) 0929 return C(y,y); 0930 0931 #if USE_CONTINUED_FRACTION 0932 ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term 0933 #else 0934 if (y < 0) { 0935 /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so 0936 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible 0937 if y*y - x*x > -36 or so. So, compute this term just in case. 0938 We also need the -exp(-x*x) term to compute Re[w] accurately 0939 in the case where y is very small. */ 0940 ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y); 0941 } 0942 else 0943 ret = exp(-x*x); // not negligible in real part if y very small 0944 #endif 0945 // (round instead of ceil as in original paper; note that x/a > 1 here) 0946 double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0 0947 double dx = a*n0 - x; 0948 sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y); 0949 sum5 = a*n0 * sum3; 0950 double exp1 = exp(4*a*dx), exp1dn = 1; 0951 int dn; 0952 for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms 0953 double np = n0 + dn, nm = n0 - dn; 0954 double tp = exp(-sqr(a*dn+dx)); 0955 double tm = tp * (exp1dn *= exp1); // trick to get tm from tp 0956 tp /= (a2*(np*np) + y*y); 0957 tm /= (a2*(nm*nm) + y*y); 0958 sum3 += tp + tm; 0959 sum5 += a * (np * tp + nm * tm); 0960 if (a * (np * tp + nm * tm) < relerr * sum5) goto finish; 0961 } 0962 while (1) { // loop over n0+dn terms only (since n0-dn <= 0) 0963 double np = n0 + dn++; 0964 double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y); 0965 sum3 += tp; 0966 sum5 += a * np * tp; 0967 if (a * np * tp < relerr * sum5) goto finish; 0968 } 0969 } 0970 finish: 0971 return ret + C((0.5*c)*y*(sum2+sum3), 0972 (0.5*c)*copysign(sum5-sum4, creal(z))); 0973 } 0974 0975 ///////////////////////////////////////////////////////////////////////// 0976 0977 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by 0978 Steven G. Johnson, October 2012. 0979 0980 This function combines a few different ideas. 0981 0982 First, for x > 50, it uses a continued-fraction expansion (same as 0983 for the Faddeeva function, but with algebraic simplifications for z=i*x). 0984 0985 Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations, 0986 but with two twists: 0987 0988 a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation, 0989 inspired by a similar transformation in the octave-forge/specfun 0990 erfcx by Soren Hauberg, results in much faster Chebyshev convergence 0991 than other simple transformations I have examined. 0992 0993 b) Instead of using a single Chebyshev polynomial for the entire 0994 [0,1] y interval, we break the interval up into 100 equal 0995 subintervals, with a switch/lookup table, and use much lower 0996 degree Chebyshev polynomials in each subinterval. This greatly 0997 improves performance in my tests. 0998 0999 For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), 1000 with the usual checks for overflow etcetera. 1001 1002 Performance-wise, it seems to be substantially faster than either 1003 the SLATEC DERFC function [or an erfcx function derived therefrom] 1004 or Cody's CALERF function (from netlib.org/specfun), while 1005 retaining near machine precision in accuracy. */ 1006 1007 /* Given y100=100*y, where y = 4/(4+x) for x >= 0, compute erfc(x). 1008 1009 Uses a look-up table of 100 different Chebyshev polynomials 1010 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated 1011 with the help of Maple and a little shell script. This allows 1012 the Chebyshev polynomials to be of significantly lower degree (about 1/4) 1013 compared to fitting the whole [0,1] interval with a single polynomial. */ 1014 static double erfcx_y100(double y100) 1015 { 1016 switch ((int) y100) { 1017 case 0: { 1018 double t = 2*y100 - 1; 1019 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; 1020 } 1021 case 1: { 1022 double t = 2*y100 - 3; 1023 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; 1024 } 1025 case 2: { 1026 double t = 2*y100 - 5; 1027 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; 1028 } 1029 case 3: { 1030 double t = 2*y100 - 7; 1031 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; 1032 } 1033 case 4: { 1034 double t = 2*y100 - 9; 1035 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; 1036 } 1037 case 5: { 1038 double t = 2*y100 - 11; 1039 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; 1040 } 1041 case 6: { 1042 double t = 2*y100 - 13; 1043 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; 1044 } 1045 case 7: { 1046 double t = 2*y100 - 15; 1047 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; 1048 } 1049 case 8: { 1050 double t = 2*y100 - 17; 1051 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; 1052 } 1053 case 9: { 1054 double t = 2*y100 - 19; 1055 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; 1056 } 1057 case 10: { 1058 double t = 2*y100 - 21; 1059 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; 1060 } 1061 case 11: { 1062 double t = 2*y100 - 23; 1063 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; 1064 } 1065 case 12: { 1066 double t = 2*y100 - 25; 1067 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; 1068 } 1069 case 13: { 1070 double t = 2*y100 - 27; 1071 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; 1072 } 1073 case 14: { 1074 double t = 2*y100 - 29; 1075 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; 1076 } 1077 case 15: { 1078 double t = 2*y100 - 31; 1079 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; 1080 } 1081 case 16: { 1082 double t = 2*y100 - 33; 1083 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; 1084 } 1085 case 17: { 1086 double t = 2*y100 - 35; 1087 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; 1088 } 1089 case 18: { 1090 double t = 2*y100 - 37; 1091 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; 1092 } 1093 case 19: { 1094 double t = 2*y100 - 39; 1095 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; 1096 } 1097 case 20: { 1098 double t = 2*y100 - 41; 1099 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; 1100 } 1101 case 21: { 1102 double t = 2*y100 - 43; 1103 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; 1104 } 1105 case 22: { 1106 double t = 2*y100 - 45; 1107 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; 1108 } 1109 case 23: { 1110 double t = 2*y100 - 47; 1111 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; 1112 } 1113 case 24: { 1114 double t = 2*y100 - 49; 1115 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; 1116 } 1117 case 25: { 1118 double t = 2*y100 - 51; 1119 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; 1120 } 1121 case 26: { 1122 double t = 2*y100 - 53; 1123 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; 1124 } 1125 case 27: { 1126 double t = 2*y100 - 55; 1127 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; 1128 } 1129 case 28: { 1130 double t = 2*y100 - 57; 1131 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; 1132 } 1133 case 29: { 1134 double t = 2*y100 - 59; 1135 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; 1136 } 1137 case 30: { 1138 double t = 2*y100 - 61; 1139 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; 1140 } 1141 case 31: { 1142 double t = 2*y100 - 63; 1143 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; 1144 } 1145 case 32: { 1146 double t = 2*y100 - 65; 1147 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; 1148 } 1149 case 33: { 1150 double t = 2*y100 - 67; 1151 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; 1152 } 1153 case 34: { 1154 double t = 2*y100 - 69; 1155 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; 1156 } 1157 case 35: { 1158 double t = 2*y100 - 71; 1159 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; 1160 } 1161 case 36: { 1162 double t = 2*y100 - 73; 1163 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; 1164 } 1165 case 37: { 1166 double t = 2*y100 - 75; 1167 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; 1168 } 1169 case 38: { 1170 double t = 2*y100 - 77; 1171 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; 1172 } 1173 case 39: { 1174 double t = 2*y100 - 79; 1175 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; 1176 } 1177 case 40: { 1178 double t = 2*y100 - 81; 1179 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; 1180 } 1181 case 41: { 1182 double t = 2*y100 - 83; 1183 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; 1184 } 1185 case 42: { 1186 double t = 2*y100 - 85; 1187 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; 1188 } 1189 case 43: { 1190 double t = 2*y100 - 87; 1191 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; 1192 } 1193 case 44: { 1194 double t = 2*y100 - 89; 1195 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; 1196 } 1197 case 45: { 1198 double t = 2*y100 - 91; 1199 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; 1200 } 1201 case 46: { 1202 double t = 2*y100 - 93; 1203 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; 1204 } 1205 case 47: { 1206 double t = 2*y100 - 95; 1207 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; 1208 } 1209 case 48: { 1210 double t = 2*y100 - 97; 1211 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; 1212 } 1213 case 49: { 1214 double t = 2*y100 - 99; 1215 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; 1216 } 1217 case 50: { 1218 double t = 2*y100 - 101; 1219 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; 1220 } 1221 case 51: { 1222 double t = 2*y100 - 103; 1223 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; 1224 } 1225 case 52: { 1226 double t = 2*y100 - 105; 1227 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; 1228 } 1229 case 53: { 1230 double t = 2*y100 - 107; 1231 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; 1232 } 1233 case 54: { 1234 double t = 2*y100 - 109; 1235 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; 1236 } 1237 case 55: { 1238 double t = 2*y100 - 111; 1239 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; 1240 } 1241 case 56: { 1242 double t = 2*y100 - 113; 1243 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; 1244 } 1245 case 57: { 1246 double t = 2*y100 - 115; 1247 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; 1248 } 1249 case 58: { 1250 double t = 2*y100 - 117; 1251 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; 1252 } 1253 case 59: { 1254 double t = 2*y100 - 119; 1255 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; 1256 } 1257 case 60: { 1258 double t = 2*y100 - 121; 1259 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; 1260 } 1261 case 61: { 1262 double t = 2*y100 - 123; 1263 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; 1264 } 1265 case 62: { 1266 double t = 2*y100 - 125; 1267 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; 1268 } 1269 case 63: { 1270 double t = 2*y100 - 127; 1271 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; 1272 } 1273 case 64: { 1274 double t = 2*y100 - 129; 1275 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; 1276 } 1277 case 65: { 1278 double t = 2*y100 - 131; 1279 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; 1280 } 1281 case 66: { 1282 double t = 2*y100 - 133; 1283 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; 1284 } 1285 case 67: { 1286 double t = 2*y100 - 135; 1287 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; 1288 } 1289 case 68: { 1290 double t = 2*y100 - 137; 1291 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; 1292 } 1293 case 69: { 1294 double t = 2*y100 - 139; 1295 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; 1296 } 1297 case 70: { 1298 double t = 2*y100 - 141; 1299 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; 1300 } 1301 case 71: { 1302 double t = 2*y100 - 143; 1303 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; 1304 } 1305 case 72: { 1306 double t = 2*y100 - 145; 1307 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; 1308 } 1309 case 73: { 1310 double t = 2*y100 - 147; 1311 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; 1312 } 1313 case 74: { 1314 double t = 2*y100 - 149; 1315 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; 1316 } 1317 case 75: { 1318 double t = 2*y100 - 151; 1319 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; 1320 } 1321 case 76: { 1322 double t = 2*y100 - 153; 1323 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; 1324 } 1325 case 77: { 1326 double t = 2*y100 - 155; 1327 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; 1328 } 1329 case 78: { 1330 double t = 2*y100 - 157; 1331 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; 1332 } 1333 case 79: { 1334 double t = 2*y100 - 159; 1335 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; 1336 } 1337 case 80: { 1338 double t = 2*y100 - 161; 1339 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; 1340 } 1341 case 81: { 1342 double t = 2*y100 - 163; 1343 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; 1344 } 1345 case 82: { 1346 double t = 2*y100 - 165; 1347 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; 1348 } 1349 case 83: { 1350 double t = 2*y100 - 167; 1351 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; 1352 } 1353 case 84: { 1354 double t = 2*y100 - 169; 1355 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; 1356 } 1357 case 85: { 1358 double t = 2*y100 - 171; 1359 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; 1360 } 1361 case 86: { 1362 double t = 2*y100 - 173; 1363 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; 1364 } 1365 case 87: { 1366 double t = 2*y100 - 175; 1367 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; 1368 } 1369 case 88: { 1370 double t = 2*y100 - 177; 1371 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; 1372 } 1373 case 89: { 1374 double t = 2*y100 - 179; 1375 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; 1376 } 1377 case 90: { 1378 double t = 2*y100 - 181; 1379 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; 1380 } 1381 case 91: { 1382 double t = 2*y100 - 183; 1383 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; 1384 } 1385 case 92: { 1386 double t = 2*y100 - 185; 1387 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; 1388 } 1389 case 93: { 1390 double t = 2*y100 - 187; 1391 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; 1392 } 1393 case 94: { 1394 double t = 2*y100 - 189; 1395 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; 1396 } 1397 case 95: { 1398 double t = 2*y100 - 191; 1399 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; 1400 } 1401 case 96: { 1402 double t = 2*y100 - 193; 1403 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; 1404 } 1405 case 97: { 1406 double t = 2*y100 - 195; 1407 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; 1408 } 1409 case 98: { 1410 double t = 2*y100 - 197; 1411 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; 1412 } 1413 case 99: { 1414 double t = 2*y100 - 199; 1415 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; 1416 } 1417 } 1418 // we only get here if y = 1, i.e. |x| < 4*eps, in which case 1419 // erfcx is within 1e-15 of 1.. 1420 return 1.0; 1421 } 1422 1423 double FADDEEVA_RE(erfcx)(double x) 1424 { 1425 if (x >= 0) { 1426 if (x > 50) { // continued-fraction expansion is faster 1427 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1428 if (x > 5e7) // 1-term expansion, important to avoid overflow 1429 return ispi / x; 1430 /* 5-term expansion (rely on compiler for CSE), simplified from: 1431 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */ 1432 return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75)); 1433 } 1434 return erfcx_y100(400/(4+x)); 1435 } 1436 else 1437 return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x) 1438 : 2*exp(x*x) - erfcx_y100(400/(4-x))); 1439 } 1440 1441 ///////////////////////////////////////////////////////////////////////// 1442 /* Compute a scaled Dawson integral 1443 FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi) 1444 equivalent to the imaginary part w(x) for real x. 1445 1446 Uses methods similar to the erfcx calculation above: continued fractions 1447 for large |x|, a lookup table of Chebyshev polynomials for smaller |x|, 1448 and finally a Taylor expansion for |x|<0.01. 1449 1450 Steven G. Johnson, October 2012. */ 1451 1452 /* Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x). 1453 1454 Uses a look-up table of 100 different Chebyshev polynomials 1455 for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated 1456 with the help of Maple and a little shell script. This allows 1457 the Chebyshev polynomials to be of significantly lower degree (about 1/30) 1458 compared to fitting the whole [0,1] interval with a single polynomial. */ 1459 static double w_im_y100(double y100, double x) { 1460 switch ((int) y100) { 1461 case 0: { 1462 double t = 2*y100 - 1; 1463 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; 1464 } 1465 case 1: { 1466 double t = 2*y100 - 3; 1467 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; 1468 } 1469 case 2: { 1470 double t = 2*y100 - 5; 1471 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; 1472 } 1473 case 3: { 1474 double t = 2*y100 - 7; 1475 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; 1476 } 1477 case 4: { 1478 double t = 2*y100 - 9; 1479 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; 1480 } 1481 case 5: { 1482 double t = 2*y100 - 11; 1483 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; 1484 } 1485 case 6: { 1486 double t = 2*y100 - 13; 1487 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; 1488 } 1489 case 7: { 1490 double t = 2*y100 - 15; 1491 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; 1492 } 1493 case 8: { 1494 double t = 2*y100 - 17; 1495 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; 1496 } 1497 case 9: { 1498 double t = 2*y100 - 19; 1499 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; 1500 } 1501 case 10: { 1502 double t = 2*y100 - 21; 1503 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; 1504 } 1505 case 11: { 1506 double t = 2*y100 - 23; 1507 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; 1508 } 1509 case 12: { 1510 double t = 2*y100 - 25; 1511 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; 1512 } 1513 case 13: { 1514 double t = 2*y100 - 27; 1515 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; 1516 } 1517 case 14: { 1518 double t = 2*y100 - 29; 1519 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; 1520 } 1521 case 15: { 1522 double t = 2*y100 - 31; 1523 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; 1524 } 1525 case 16: { 1526 double t = 2*y100 - 33; 1527 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; 1528 } 1529 case 17: { 1530 double t = 2*y100 - 35; 1531 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; 1532 } 1533 case 18: { 1534 double t = 2*y100 - 37; 1535 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; 1536 } 1537 case 19: { 1538 double t = 2*y100 - 39; 1539 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; 1540 } 1541 case 20: { 1542 double t = 2*y100 - 41; 1543 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; 1544 } 1545 case 21: { 1546 double t = 2*y100 - 43; 1547 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; 1548 } 1549 case 22: { 1550 double t = 2*y100 - 45; 1551 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; 1552 } 1553 case 23: { 1554 double t = 2*y100 - 47; 1555 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; 1556 } 1557 case 24: { 1558 double t = 2*y100 - 49; 1559 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; 1560 } 1561 case 25: { 1562 double t = 2*y100 - 51; 1563 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; 1564 } 1565 case 26: { 1566 double t = 2*y100 - 53; 1567 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; 1568 } 1569 case 27: { 1570 double t = 2*y100 - 55; 1571 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; 1572 } 1573 case 28: { 1574 double t = 2*y100 - 57; 1575 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; 1576 } 1577 case 29: { 1578 double t = 2*y100 - 59; 1579 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; 1580 } 1581 case 30: { 1582 double t = 2*y100 - 61; 1583 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; 1584 } 1585 case 31: { 1586 double t = 2*y100 - 63; 1587 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; 1588 } 1589 case 32: { 1590 double t = 2*y100 - 65; 1591 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; 1592 } 1593 case 33: { 1594 double t = 2*y100 - 67; 1595 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; 1596 } 1597 case 34: { 1598 double t = 2*y100 - 69; 1599 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; 1600 } 1601 case 35: { 1602 double t = 2*y100 - 71; 1603 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; 1604 } 1605 case 36: { 1606 double t = 2*y100 - 73; 1607 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; 1608 } 1609 case 37: { 1610 double t = 2*y100 - 75; 1611 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; 1612 } 1613 case 38: { 1614 double t = 2*y100 - 77; 1615 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; 1616 } 1617 case 39: { 1618 double t = 2*y100 - 79; 1619 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; 1620 } 1621 case 40: { 1622 double t = 2*y100 - 81; 1623 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; 1624 } 1625 case 41: { 1626 double t = 2*y100 - 83; 1627 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; 1628 } 1629 case 42: { 1630 double t = 2*y100 - 85; 1631 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; 1632 } 1633 case 43: { 1634 double t = 2*y100 - 87; 1635 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; 1636 } 1637 case 44: { 1638 double t = 2*y100 - 89; 1639 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; 1640 } 1641 case 45: { 1642 double t = 2*y100 - 91; 1643 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; 1644 } 1645 case 46: { 1646 double t = 2*y100 - 93; 1647 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; 1648 } 1649 case 47: { 1650 double t = 2*y100 - 95; 1651 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; 1652 } 1653 case 48: { 1654 double t = 2*y100 - 97; 1655 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; 1656 } 1657 case 49: { 1658 double t = 2*y100 - 99; 1659 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; 1660 } 1661 case 50: { 1662 double t = 2*y100 - 101; 1663 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; 1664 } 1665 case 51: { 1666 double t = 2*y100 - 103; 1667 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; 1668 } 1669 case 52: { 1670 double t = 2*y100 - 105; 1671 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; 1672 } 1673 case 53: { 1674 double t = 2*y100 - 107; 1675 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; 1676 } 1677 case 54: { 1678 double t = 2*y100 - 109; 1679 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; 1680 } 1681 case 55: { 1682 double t = 2*y100 - 111; 1683 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; 1684 } 1685 case 56: { 1686 double t = 2*y100 - 113; 1687 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; 1688 } 1689 case 57: { 1690 double t = 2*y100 - 115; 1691 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; 1692 } 1693 case 58: { 1694 double t = 2*y100 - 117; 1695 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; 1696 } 1697 case 59: { 1698 double t = 2*y100 - 119; 1699 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; 1700 } 1701 case 60: { 1702 double t = 2*y100 - 121; 1703 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; 1704 } 1705 case 61: { 1706 double t = 2*y100 - 123; 1707 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; 1708 } 1709 case 62: { 1710 double t = 2*y100 - 125; 1711 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; 1712 } 1713 case 63: { 1714 double t = 2*y100 - 127; 1715 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; 1716 } 1717 case 64: { 1718 double t = 2*y100 - 129; 1719 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; 1720 } 1721 case 65: { 1722 double t = 2*y100 - 131; 1723 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; 1724 } 1725 case 66: { 1726 double t = 2*y100 - 133; 1727 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; 1728 } 1729 case 67: { 1730 double t = 2*y100 - 135; 1731 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; 1732 } 1733 case 68: { 1734 double t = 2*y100 - 137; 1735 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; 1736 } 1737 case 69: { 1738 double t = 2*y100 - 139; 1739 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; 1740 } 1741 case 70: { 1742 double t = 2*y100 - 141; 1743 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; 1744 } 1745 case 71: { 1746 double t = 2*y100 - 143; 1747 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; 1748 } 1749 case 72: { 1750 double t = 2*y100 - 145; 1751 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; 1752 } 1753 case 73: { 1754 double t = 2*y100 - 147; 1755 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; 1756 } 1757 case 74: { 1758 double t = 2*y100 - 149; 1759 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; 1760 } 1761 case 75: { 1762 double t = 2*y100 - 151; 1763 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; 1764 } 1765 case 76: { 1766 double t = 2*y100 - 153; 1767 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; 1768 } 1769 case 77: { 1770 double t = 2*y100 - 155; 1771 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; 1772 } 1773 case 78: { 1774 double t = 2*y100 - 157; 1775 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; 1776 } 1777 case 79: { 1778 double t = 2*y100 - 159; 1779 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; 1780 } 1781 case 80: { 1782 double t = 2*y100 - 161; 1783 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; 1784 } 1785 case 81: { 1786 double t = 2*y100 - 163; 1787 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; 1788 } 1789 case 82: { 1790 double t = 2*y100 - 165; 1791 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; 1792 } 1793 case 83: { 1794 double t = 2*y100 - 167; 1795 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; 1796 } 1797 case 84: { 1798 double t = 2*y100 - 169; 1799 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; 1800 } 1801 case 85: { 1802 double t = 2*y100 - 171; 1803 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; 1804 } 1805 case 86: { 1806 double t = 2*y100 - 173; 1807 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; 1808 } 1809 case 87: { 1810 double t = 2*y100 - 175; 1811 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; 1812 } 1813 case 88: { 1814 double t = 2*y100 - 177; 1815 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; 1816 } 1817 case 89: { 1818 double t = 2*y100 - 179; 1819 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; 1820 } 1821 case 90: { 1822 double t = 2*y100 - 181; 1823 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; 1824 } 1825 case 91: { 1826 double t = 2*y100 - 183; 1827 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; 1828 } 1829 case 92: { 1830 double t = 2*y100 - 185; 1831 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; 1832 } 1833 case 93: { 1834 double t = 2*y100 - 187; 1835 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; 1836 } 1837 case 94: { 1838 double t = 2*y100 - 189; 1839 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; 1840 } 1841 case 95: { 1842 double t = 2*y100 - 191; 1843 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; 1844 } 1845 case 96: { 1846 double t = 2*y100 - 193; 1847 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; 1848 } 1849 case 97: case 98: 1850 case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...) 1851 // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9) 1852 double x2 = x*x; 1853 return x * (1.1283791670955125739 1854 - x2 * (0.75225277806367504925 1855 - x2 * (0.30090111122547001970 1856 - x2 * (0.085971746064420005629 1857 - x2 * 0.016931216931216931217)))); 1858 } 1859 } 1860 /* Since 0 <= y100 < 101, this is only reached if x is NaN, 1861 in which case we should return NaN. */ 1862 return NaN; 1863 } 1864 1865 double FADDEEVA(w_im)(double x) 1866 { 1867 if (x >= 0) { 1868 if (x > 45) { // continued-fraction expansion is faster 1869 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1870 if (x > 5e7) // 1-term expansion, important to avoid overflow 1871 return ispi / x; 1872 /* 5-term expansion (rely on compiler for CSE), simplified from: 1873 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ 1874 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75)); 1875 } 1876 return w_im_y100(100/(1+x), x); 1877 } 1878 else { // = -FADDEEVA(w_im)(-x) 1879 if (x < -45) { // continued-fraction expansion is faster 1880 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1881 if (x < -5e7) // 1-term expansion, important to avoid overflow 1882 return ispi / x; 1883 /* 5-term expansion (rely on compiler for CSE), simplified from: 1884 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ 1885 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75)); 1886 } 1887 return -w_im_y100(100/(1-x), -x); 1888 } 1889 } 1890 1891 ///////////////////////////////////////////////////////////////////////// 1892 1893 // Compile with -DTEST_FADDEEVA to compile a little test program 1894 #ifdef TEST_FADDEEVA 1895 1896 #ifdef __cplusplus 1897 # include <cstdio> 1898 #else 1899 # include <stdio.h> 1900 #endif 1901 1902 // compute relative error |b-a|/|a|, handling case of NaN and Inf, 1903 static double relerr(double a, double b) { 1904 if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) { 1905 if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) || 1906 (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) || 1907 (isinf(a) && isinf(b) && a*b < 0)) 1908 return Inf; // "infinite" error 1909 return 0; // matching infinity/nan results counted as zero error 1910 } 1911 if (a == 0) 1912 return b == 0 ? 0 : Inf; 1913 else 1914 return fabs((b-a) / a); 1915 } 1916 1917 int main(void) { 1918 double errmax_all = 0; 1919 { 1920 printf("############# w(z) tests #############\n"); 1921 #define NTST 57 // define instead of const for C compatibility 1922 cmplx z[NTST] = { 1923 C(624.2,-0.26123), 1924 C(-0.4,3.), 1925 C(0.6,2.), 1926 C(-1.,1.), 1927 C(-1.,-9.), 1928 C(-1.,9.), 1929 C(-0.0000000234545,1.1234), 1930 C(-3.,5.1), 1931 C(-53,30.1), 1932 C(0.0,0.12345), 1933 C(11,1), 1934 C(-22,-2), 1935 C(9,-28), 1936 C(21,-33), 1937 C(1e5,1e5), 1938 C(1e14,1e14), 1939 C(-3001,-1000), 1940 C(1e160,-1e159), 1941 C(-6.01,0.01), 1942 C(-0.7,-0.7), 1943 C(2.611780000000000e+01, 4.540909610972489e+03), 1944 C(0.8e7,0.3e7), 1945 C(-20,-19.8081), 1946 C(1e-16,-1.1e-16), 1947 C(2.3e-8,1.3e-8), 1948 C(6.3,-1e-13), 1949 C(6.3,1e-20), 1950 C(1e-20,6.3), 1951 C(1e-20,16.3), 1952 C(9,1e-300), 1953 C(6.01,0.11), 1954 C(8.01,1.01e-10), 1955 C(28.01,1e-300), 1956 C(10.01,1e-200), 1957 C(10.01,-1e-200), 1958 C(10.01,0.99e-10), 1959 C(10.01,-0.99e-10), 1960 C(1e-20,7.01), 1961 C(-1,7.01), 1962 C(5.99,7.01), 1963 C(1,0), 1964 C(55,0), 1965 C(-0.1,0), 1966 C(1e-20,0), 1967 C(0,5e-14), 1968 C(0,51), 1969 C(Inf,0), 1970 C(-Inf,0), 1971 C(0,Inf), 1972 C(0,-Inf), 1973 C(Inf,Inf), 1974 C(Inf,-Inf), 1975 C(NaN,NaN), 1976 C(NaN,0), 1977 C(0,NaN), 1978 C(NaN,Inf), 1979 C(Inf,NaN) 1980 }; 1981 cmplx w[NTST] = { /* w(z), computed with WolframAlpha 1982 ... note that WolframAlpha is problematic 1983 some of the above inputs, so I had to 1984 use the continued-fraction expansion 1985 in WolframAlpha in some cases, or switch 1986 to Maple */ 1987 C(-3.78270245518980507452677445620103199303131110e-7, 1988 0.000903861276433172057331093754199933411710053155), 1989 C(0.1764906227004816847297495349730234591778719532788, 1990 -0.02146550539468457616788719893991501311573031095617), 1991 C(0.2410250715772692146133539023007113781272362309451, 1992 0.06087579663428089745895459735240964093522265589350), 1993 C(0.30474420525691259245713884106959496013413834051768, 1994 -0.20821893820283162728743734725471561394145872072738), 1995 C(7.317131068972378096865595229600561710140617977e34, 1996 8.321873499714402777186848353320412813066170427e34), 1997 C(0.0615698507236323685519612934241429530190806818395, 1998 -0.00676005783716575013073036218018565206070072304635), 1999 C(0.3960793007699874918961319170187598400134746631, 2000 -5.593152259116644920546186222529802777409274656e-9), 2001 C(0.08217199226739447943295069917990417630675021771804, 2002 -0.04701291087643609891018366143118110965272615832184), 2003 C(0.00457246000350281640952328010227885008541748668738, 2004 -0.00804900791411691821818731763401840373998654987934), 2005 C(0.8746342859608052666092782112565360755791467973338452, 2006 0.), 2007 C(0.00468190164965444174367477874864366058339647648741, 2008 0.0510735563901306197993676329845149741675029197050), 2009 C(-0.0023193175200187620902125853834909543869428763219, 2010 -0.025460054739731556004902057663500272721780776336), 2011 C(9.11463368405637174660562096516414499772662584e304, 2012 3.97101807145263333769664875189354358563218932e305), 2013 C(-4.4927207857715598976165541011143706155432296e281, 2014 -2.8019591213423077494444700357168707775769028e281), 2015 C(2.820947917809305132678577516325951485807107151e-6, 2016 2.820947917668257736791638444590253942253354058e-6), 2017 C(2.82094791773878143474039725787438662716372268e-15, 2018 2.82094791773878143474039725773333923127678361e-15), 2019 C(-0.0000563851289696244350147899376081488003110150498, 2020 -0.000169211755126812174631861529808288295454992688), 2021 C(-5.586035480670854326218608431294778077663867e-162, 2022 5.586035480670854326218608431294778077663867e-161), 2023 C(0.00016318325137140451888255634399123461580248456, 2024 -0.095232456573009287370728788146686162555021209999), 2025 C(0.69504753678406939989115375989939096800793577783885, 2026 -1.8916411171103639136680830887017670616339912024317), 2027 C(0.0001242418269653279656612334210746733213167234822, 2028 7.145975826320186888508563111992099992116786763e-7), 2029 C(2.318587329648353318615800865959225429377529825e-8, 2030 6.182899545728857485721417893323317843200933380e-8), 2031 C(-0.0133426877243506022053521927604277115767311800303, 2032 -0.0148087097143220769493341484176979826888871576145), 2033 C(1.00000000000000012412170838050638522857747934, 2034 1.12837916709551279389615890312156495593616433e-16), 2035 C(0.9999999853310704677583504063775310832036830015, 2036 2.595272024519678881897196435157270184030360773e-8), 2037 C(-1.4731421795638279504242963027196663601154624e-15, 2038 0.090727659684127365236479098488823462473074709), 2039 C(5.79246077884410284575834156425396800754409308e-18, 2040 0.0907276596841273652364790985059772809093822374), 2041 C(0.0884658993528521953466533278764830881245144368, 2042 1.37088352495749125283269718778582613192166760e-22), 2043 C(0.0345480845419190424370085249304184266813447878, 2044 2.11161102895179044968099038990446187626075258e-23), 2045 C(6.63967719958073440070225527042829242391918213e-36, 2046 0.0630820900592582863713653132559743161572639353), 2047 C(0.00179435233208702644891092397579091030658500743634, 2048 0.0951983814805270647939647438459699953990788064762), 2049 C(9.09760377102097999924241322094863528771095448e-13, 2050 0.0709979210725138550986782242355007611074966717), 2051 C(7.2049510279742166460047102593255688682910274423e-304, 2052 0.0201552956479526953866611812593266285000876784321), 2053 C(3.04543604652250734193622967873276113872279682e-44, 2054 0.0566481651760675042930042117726713294607499165), 2055 C(3.04543604652250734193622967873276113872279682e-44, 2056 0.0566481651760675042930042117726713294607499165), 2057 C(0.5659928732065273429286988428080855057102069081e-12, 2058 0.056648165176067504292998527162143030538756683302), 2059 C(-0.56599287320652734292869884280802459698927645e-12, 2060 0.0566481651760675042929985271621430305387566833029), 2061 C(0.0796884251721652215687859778119964009569455462, 2062 1.11474461817561675017794941973556302717225126e-22), 2063 C(0.07817195821247357458545539935996687005781943386550, 2064 -0.01093913670103576690766705513142246633056714279654), 2065 C(0.04670032980990449912809326141164730850466208439937, 2066 0.03944038961933534137558064191650437353429669886545), 2067 C(0.36787944117144232159552377016146086744581113103176, 2068 0.60715770584139372911503823580074492116122092866515), 2069 C(0, 2070 0.010259688805536830986089913987516716056946786526145), 2071 C(0.99004983374916805357390597718003655777207908125383, 2072 -0.11208866436449538036721343053869621153527769495574), 2073 C(0.99999999999999999999999999999999999999990000, 2074 1.12837916709551257389615890312154517168802603e-20), 2075 C(0.999999999999943581041645226871305192054749891144158, 2076 0), 2077 C(0.0110604154853277201542582159216317923453996211744250, 2078 0), 2079 C(0,0), 2080 C(0,0), 2081 C(0,0), 2082 C(Inf,0), 2083 C(0,0), 2084 C(NaN,NaN), 2085 C(NaN,NaN), 2086 C(NaN,NaN), 2087 C(NaN,0), 2088 C(NaN,NaN), 2089 C(NaN,NaN) 2090 }; 2091 double errmax = 0; 2092 for (int i = 0; i < NTST; ++i) { 2093 cmplx fw = FADDEEVA(w)(z[i],0.); 2094 double re_err = relerr(creal(w[i]), creal(fw)); 2095 double im_err = relerr(cimag(w[i]), cimag(fw)); 2096 printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", 2097 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), 2098 re_err, im_err); 2099 if (re_err > errmax) errmax = re_err; 2100 if (im_err > errmax) errmax = im_err; 2101 } 2102 if (errmax > 1e-13) { 2103 printf("FAILURE -- relative error %g too large!\n", errmax); 2104 return 1; 2105 } 2106 printf("SUCCESS (max relative error = %g)\n", errmax); 2107 if (errmax > errmax_all) errmax_all = errmax; 2108 } 2109 { 2110 #undef NTST 2111 #define NTST 41 // define instead of const for C compatibility 2112 cmplx z[NTST] = { 2113 C(1,2), 2114 C(-1,2), 2115 C(1,-2), 2116 C(-1,-2), 2117 C(9,-28), 2118 C(21,-33), 2119 C(1e3,1e3), 2120 C(-3001,-1000), 2121 C(1e160,-1e159), 2122 C(5.1e-3, 1e-8), 2123 C(-4.9e-3, 4.95e-3), 2124 C(4.9e-3, 0.5), 2125 C(4.9e-4, -0.5e1), 2126 C(-4.9e-5, -0.5e2), 2127 C(5.1e-3, 0.5), 2128 C(5.1e-4, -0.5e1), 2129 C(-5.1e-5, -0.5e2), 2130 C(1e-6,2e-6), 2131 C(0,2e-6), 2132 C(0,2), 2133 C(0,20), 2134 C(0,200), 2135 C(Inf,0), 2136 C(-Inf,0), 2137 C(0,Inf), 2138 C(0,-Inf), 2139 C(Inf,Inf), 2140 C(Inf,-Inf), 2141 C(NaN,NaN), 2142 C(NaN,0), 2143 C(0,NaN), 2144 C(NaN,Inf), 2145 C(Inf,NaN), 2146 C(1e-3,NaN), 2147 C(7e-2,7e-2), 2148 C(7e-2,-7e-4), 2149 C(-9e-2,7e-4), 2150 C(-9e-2,9e-2), 2151 C(-7e-4,9e-2), 2152 C(7e-2,0.9e-2), 2153 C(7e-2,1.1e-2) 2154 }; 2155 cmplx w[NTST] = { // erf(z[i]), evaluated with Maple 2156 C(-0.5366435657785650339917955593141927494421, 2157 -5.049143703447034669543036958614140565553), 2158 C(0.5366435657785650339917955593141927494421, 2159 -5.049143703447034669543036958614140565553), 2160 C(-0.5366435657785650339917955593141927494421, 2161 5.049143703447034669543036958614140565553), 2162 C(0.5366435657785650339917955593141927494421, 2163 5.049143703447034669543036958614140565553), 2164 C(0.3359473673830576996788000505817956637777e304, 2165 -0.1999896139679880888755589794455069208455e304), 2166 C(0.3584459971462946066523939204836760283645e278, 2167 0.3818954885257184373734213077678011282505e280), 2168 C(0.9996020422657148639102150147542224526887, 2169 0.00002801044116908227889681753993542916894856), 2170 C(-1, 0), 2171 C(1, 0), 2172 C(0.005754683859034800134412990541076554934877, 2173 0.1128349818335058741511924929801267822634e-7), 2174 C(-0.005529149142341821193633460286828381876955, 2175 0.005585388387864706679609092447916333443570), 2176 C(0.007099365669981359632319829148438283865814, 2177 0.6149347012854211635026981277569074001219), 2178 C(0.3981176338702323417718189922039863062440e8, 2179 -0.8298176341665249121085423917575122140650e10), 2180 C(-Inf, 2181 -Inf), 2182 C(0.007389128308257135427153919483147229573895, 2183 0.6149332524601658796226417164791221815139), 2184 C(0.4143671923267934479245651547534414976991e8, 2185 -0.8298168216818314211557046346850921446950e10), 2186 C(-Inf, 2187 -Inf), 2188 C(0.1128379167099649964175513742247082845155e-5, 2189 0.2256758334191777400570377193451519478895e-5), 2190 C(0, 2191 0.2256758334194034158904576117253481476197e-5), 2192 C(0, 2193 18.56480241457555259870429191324101719886), 2194 C(0, 2195 0.1474797539628786202447733153131835124599e173), 2196 C(0, 2197 Inf), 2198 C(1,0), 2199 C(-1,0), 2200 C(0,Inf), 2201 C(0,-Inf), 2202 C(NaN,NaN), 2203 C(NaN,NaN), 2204 C(NaN,NaN), 2205 C(NaN,0), 2206 C(0,NaN), 2207 C(NaN,NaN), 2208 C(NaN,NaN), 2209 C(NaN,NaN), 2210 C(0.07924380404615782687930591956705225541145, 2211 0.07872776218046681145537914954027729115247), 2212 C(0.07885775828512276968931773651224684454495, 2213 -0.0007860046704118224342390725280161272277506), 2214 C(-0.1012806432747198859687963080684978759881, 2215 0.0007834934747022035607566216654982820299469), 2216 C(-0.1020998418798097910247132140051062512527, 2217 0.1010030778892310851309082083238896270340), 2218 C(-0.0007962891763147907785684591823889484764272, 2219 0.1018289385936278171741809237435404896152), 2220 C(0.07886408666470478681566329888615410479530, 2221 0.01010604288780868961492224347707949372245), 2222 C(0.07886723099940260286824654364807981336591, 2223 0.01235199327873258197931147306290916629654) 2224 }; 2225 #define TST(f,isc) \ 2226 printf("############# " #f "(z) tests #############\n"); \ 2227 double errmax = 0; \ 2228 for (int i = 0; i < NTST; ++i) { \ 2229 cmplx fw = FADDEEVA(f)(z[i],0.); \ 2230 double re_err = relerr(creal(w[i]), creal(fw)); \ 2231 double im_err = relerr(cimag(w[i]), cimag(fw)); \ 2232 printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \ 2233 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \ 2234 re_err, im_err); \ 2235 if (re_err > errmax) errmax = re_err; \ 2236 if (im_err > errmax) errmax = im_err; \ 2237 } \ 2238 if (errmax > 1e-13) { \ 2239 printf("FAILURE -- relative error %g too large!\n", errmax); \ 2240 return 1; \ 2241 } \ 2242 printf("Checking " #f "(x) special case...\n"); \ 2243 for (int i = 0; i < 10000; ++i) { \ 2244 double x = pow(10., -300. + i * 600. / (10000 - 1)); \ 2245 double re_err = relerr(FADDEEVA_RE(f)(x), \ 2246 creal(FADDEEVA(f)(C(x,x*isc),0.))); \ 2247 if (re_err > errmax) errmax = re_err; \ 2248 re_err = relerr(FADDEEVA_RE(f)(-x), \ 2249 creal(FADDEEVA(f)(C(-x,x*isc),0.))); \ 2250 if (re_err > errmax) errmax = re_err; \ 2251 } \ 2252 { \ 2253 double re_err = relerr(FADDEEVA_RE(f)(Inf), \ 2254 creal(FADDEEVA(f)(C(Inf,0.),0.))); \ 2255 if (re_err > errmax) errmax = re_err; \ 2256 re_err = relerr(FADDEEVA_RE(f)(-Inf), \ 2257 creal(FADDEEVA(f)(C(-Inf,0.),0.))); \ 2258 if (re_err > errmax) errmax = re_err; \ 2259 re_err = relerr(FADDEEVA_RE(f)(NaN), \ 2260 creal(FADDEEVA(f)(C(NaN,0.),0.))); \ 2261 if (re_err > errmax) errmax = re_err; \ 2262 } \ 2263 if (errmax > 1e-13) { \ 2264 printf("FAILURE -- relative error %g too large!\n", errmax); \ 2265 return 1; \ 2266 } \ 2267 printf("SUCCESS (max relative error = %g)\n", errmax); \ 2268 if (errmax > errmax_all) errmax_all = errmax 2269 2270 TST(erf, 1e-20); 2271 } 2272 { 2273 // since erfi just calls through to erf, just one test should 2274 // be sufficient to make sure I didn't screw up the signs or something 2275 #undef NTST 2276 #define NTST 1 // define instead of const for C compatibility 2277 cmplx z[NTST] = { C(1.234,0.5678) }; 2278 cmplx w[NTST] = { // erfi(z[i]), computed with Maple 2279 C(1.081032284405373149432716643834106923212, 2280 1.926775520840916645838949402886591180834) 2281 }; 2282 TST(erfi, 0); 2283 } 2284 { 2285 // since erfcx just calls through to w, just one test should 2286 // be sufficient to make sure I didn't screw up the signs or something 2287 #undef NTST 2288 #define NTST 1 // define instead of const for C compatibility 2289 cmplx z[NTST] = { C(1.234,0.5678) }; 2290 cmplx w[NTST] = { // erfcx(z[i]), computed with Maple 2291 C(0.3382187479799972294747793561190487832579, 2292 -0.1116077470811648467464927471872945833154) 2293 }; 2294 TST(erfcx, 0); 2295 } 2296 { 2297 #undef NTST 2298 #define NTST 30 // define instead of const for C compatibility 2299 cmplx z[NTST] = { 2300 C(1,2), 2301 C(-1,2), 2302 C(1,-2), 2303 C(-1,-2), 2304 C(9,-28), 2305 C(21,-33), 2306 C(1e3,1e3), 2307 C(-3001,-1000), 2308 C(1e160,-1e159), 2309 C(5.1e-3, 1e-8), 2310 C(0,2e-6), 2311 C(0,2), 2312 C(0,20), 2313 C(0,200), 2314 C(2e-6,0), 2315 C(2,0), 2316 C(20,0), 2317 C(200,0), 2318 C(Inf,0), 2319 C(-Inf,0), 2320 C(0,Inf), 2321 C(0,-Inf), 2322 C(Inf,Inf), 2323 C(Inf,-Inf), 2324 C(NaN,NaN), 2325 C(NaN,0), 2326 C(0,NaN), 2327 C(NaN,Inf), 2328 C(Inf,NaN), 2329 C(88,0) 2330 }; 2331 cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple 2332 C(1.536643565778565033991795559314192749442, 2333 5.049143703447034669543036958614140565553), 2334 C(0.4633564342214349660082044406858072505579, 2335 5.049143703447034669543036958614140565553), 2336 C(1.536643565778565033991795559314192749442, 2337 -5.049143703447034669543036958614140565553), 2338 C(0.4633564342214349660082044406858072505579, 2339 -5.049143703447034669543036958614140565553), 2340 C(-0.3359473673830576996788000505817956637777e304, 2341 0.1999896139679880888755589794455069208455e304), 2342 C(-0.3584459971462946066523939204836760283645e278, 2343 -0.3818954885257184373734213077678011282505e280), 2344 C(0.0003979577342851360897849852457775473112748, 2345 -0.00002801044116908227889681753993542916894856), 2346 C(2, 0), 2347 C(0, 0), 2348 C(0.9942453161409651998655870094589234450651, 2349 -0.1128349818335058741511924929801267822634e-7), 2350 C(1, 2351 -0.2256758334194034158904576117253481476197e-5), 2352 C(1, 2353 -18.56480241457555259870429191324101719886), 2354 C(1, 2355 -0.1474797539628786202447733153131835124599e173), 2356 C(1, -Inf), 2357 C(0.9999977432416658119838633199332831406314, 2358 0), 2359 C(0.004677734981047265837930743632747071389108, 2360 0), 2361 C(0.5395865611607900928934999167905345604088e-175, 2362 0), 2363 C(0, 0), 2364 C(0, 0), 2365 C(2, 0), 2366 C(1, -Inf), 2367 C(1, Inf), 2368 C(NaN, NaN), 2369 C(NaN, NaN), 2370 C(NaN, NaN), 2371 C(NaN, 0), 2372 C(1, NaN), 2373 C(NaN, NaN), 2374 C(NaN, NaN), 2375 C(0,0) 2376 }; 2377 TST(erfc, 1e-20); 2378 } 2379 { 2380 #undef NTST 2381 #define NTST 48 // define instead of const for C compatibility 2382 cmplx z[NTST] = { 2383 C(2,1), 2384 C(-2,1), 2385 C(2,-1), 2386 C(-2,-1), 2387 C(-28,9), 2388 C(33,-21), 2389 C(1e3,1e3), 2390 C(-1000,-3001), 2391 C(1e-8, 5.1e-3), 2392 C(4.95e-3, -4.9e-3), 2393 C(5.1e-3, 5.1e-3), 2394 C(0.5, 4.9e-3), 2395 C(-0.5e1, 4.9e-4), 2396 C(-0.5e2, -4.9e-5), 2397 C(0.5e3, 4.9e-6), 2398 C(0.5, 5.1e-3), 2399 C(-0.5e1, 5.1e-4), 2400 C(-0.5e2, -5.1e-5), 2401 C(1e-6,2e-6), 2402 C(2e-6,0), 2403 C(2,0), 2404 C(20,0), 2405 C(200,0), 2406 C(0,4.9e-3), 2407 C(0,-5.1e-3), 2408 C(0,2e-6), 2409 C(0,-2), 2410 C(0,20), 2411 C(0,-200), 2412 C(Inf,0), 2413 C(-Inf,0), 2414 C(0,Inf), 2415 C(0,-Inf), 2416 C(Inf,Inf), 2417 C(Inf,-Inf), 2418 C(NaN,NaN), 2419 C(NaN,0), 2420 C(0,NaN), 2421 C(NaN,Inf), 2422 C(Inf,NaN), 2423 C(39, 6.4e-5), 2424 C(41, 6.09e-5), 2425 C(4.9e7, 5e-11), 2426 C(5.1e7, 4.8e-11), 2427 C(1e9, 2.4e-12), 2428 C(1e11, 2.4e-14), 2429 C(1e13, 2.4e-16), 2430 C(1e300, 2.4e-303) 2431 }; 2432 cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple 2433 C(0.1635394094345355614904345232875688576839, 2434 -0.1531245755371229803585918112683241066853), 2435 C(-0.1635394094345355614904345232875688576839, 2436 -0.1531245755371229803585918112683241066853), 2437 C(0.1635394094345355614904345232875688576839, 2438 0.1531245755371229803585918112683241066853), 2439 C(-0.1635394094345355614904345232875688576839, 2440 0.1531245755371229803585918112683241066853), 2441 C(-0.01619082256681596362895875232699626384420, 2442 -0.005210224203359059109181555401330902819419), 2443 C(0.01078377080978103125464543240346760257008, 2444 0.006866888783433775382193630944275682670599), 2445 C(-0.5808616819196736225612296471081337245459, 2446 0.6688593905505562263387760667171706325749), 2447 C(Inf, 2448 -Inf), 2449 C(0.1000052020902036118082966385855563526705e-7, 2450 0.005100088434920073153418834680320146441685), 2451 C(0.004950156837581592745389973960217444687524, 2452 -0.004899838305155226382584756154100963570500), 2453 C(0.005100176864319675957314822982399286703798, 2454 0.005099823128319785355949825238269336481254), 2455 C(0.4244534840871830045021143490355372016428, 2456 0.002820278933186814021399602648373095266538), 2457 C(-0.1021340733271046543881236523269967674156, 2458 -0.00001045696456072005761498961861088944159916), 2459 C(-0.01000200120119206748855061636187197886859, 2460 0.9805885888237419500266621041508714123763e-8), 2461 C(0.001000002000012000023960527532953151819595, 2462 -0.9800058800588007290937355024646722133204e-11), 2463 C(0.4244549085628511778373438768121222815752, 2464 0.002935393851311701428647152230552122898291), 2465 C(-0.1021340732357117208743299813648493928105, 2466 -0.00001088377943049851799938998805451564893540), 2467 C(-0.01000200120119126652710792390331206563616, 2468 0.1020612612857282306892368985525393707486e-7), 2469 C(0.1000000000007333333333344266666666664457e-5, 2470 0.2000000000001333333333323199999999978819e-5), 2471 C(0.1999999999994666666666675199999999990248e-5, 2472 0), 2473 C(0.3013403889237919660346644392864226952119, 2474 0), 2475 C(0.02503136792640367194699495234782353186858, 2476 0), 2477 C(0.002500031251171948248596912483183760683918, 2478 0), 2479 C(0,0.004900078433419939164774792850907128053308), 2480 C(0,-0.005100088434920074173454208832365950009419), 2481 C(0,0.2000000000005333333333341866666666676419e-5), 2482 C(0,-48.16001211429122974789822893525016528191), 2483 C(0,0.4627407029504443513654142715903005954668e174), 2484 C(0,-Inf), 2485 C(0,0), 2486 C(-0,0), 2487 C(0, Inf), 2488 C(0, -Inf), 2489 C(NaN, NaN), 2490 C(NaN, NaN), 2491 C(NaN, NaN), 2492 C(NaN, 0), 2493 C(0, NaN), 2494 C(NaN, NaN), 2495 C(NaN, NaN), 2496 C(0.01282473148489433743567240624939698290584, 2497 -0.2105957276516618621447832572909153498104e-7), 2498 C(0.01219875253423634378984109995893708152885, 2499 -0.1813040560401824664088425926165834355953e-7), 2500 C(0.1020408163265306334945473399689037886997e-7, 2501 -0.1041232819658476285651490827866174985330e-25), 2502 C(0.9803921568627452865036825956835185367356e-8, 2503 -0.9227220299884665067601095648451913375754e-26), 2504 C(0.5000000000000000002500000000000000003750e-9, 2505 -0.1200000000000000001800000188712838420241e-29), 2506 C(5.00000000000000000000025000000000000000000003e-12, 2507 -1.20000000000000000000018000000000000000000004e-36), 2508 C(5.00000000000000000000000002500000000000000000e-14, 2509 -1.20000000000000000000000001800000000000000000e-42), 2510 C(5e-301, 0) 2511 }; 2512 TST(Dawson, 1e-20); 2513 } 2514 printf("#####################################\n"); 2515 printf("SUCCESS (max relative error = %g)\n", errmax_all); 2516 } 2517 2518 #endif