Mercurial > octave-libgccjit

comparison liboctave/cruft/Faddeeva/Faddeeva.cc @ 15801:2e30a1aadff2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge new upstream Faddeeva release. * Faddeeva.cc: Use numeric_limits for Inf and NaN, support C11 isnan/isinf, use gnulib floor and copysign, convert tabs to spaces, slight accuracy improvements to erf and dawson functions for |z| around 1e-2.
author Steven G. Johnson <stevenj@alum.mit.edu>
date Sat, 15 Dec 2012 14:12:47 -0500
parents 2fac72a256ce
children 161ebb78ac1b
comparison
equal deleted inserted replaced
15793:06832c90ae7d 15801:2e30a1aadff2
1 // -*- mode:c++; tab-width:2; indent-tabs-mode:nil; -*-
2
1 /* Copyright (c) 2012 Massachusetts Institute of Technology 3 /* Copyright (c) 2012 Massachusetts Institute of Technology
2 * 4 *
3 * Permission is hereby granted, free of charge, to any person obtaining 5 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the 6 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including 7 * "Software"), to deal in the Software without restriction, including
96 98
97 REVISION HISTORY: 99 REVISION HISTORY:
98 4 October 2012: Initial public release (SGJ) 100 4 October 2012: Initial public release (SGJ)
99 5 October 2012: Revised (SGJ) to fix spelling error, 101 5 October 2012: Revised (SGJ) to fix spelling error,
100 start summation for large x at round(x/a) (> 1) 102 start summation for large x at round(x/a) (> 1)
101 rather than ceil(x/a) as in the original 103 rather than ceil(x/a) as in the original
102 paper, which should slightly improve performance 104 paper, which should slightly improve performance
103 (and, apparently, slightly improves accuracy) 105 (and, apparently, slightly improves accuracy)
104 19 October 2012: Revised (SGJ) to fix bugs for large x, large -y, 106 19 October 2012: Revised (SGJ) to fix bugs for large x, large -y,
105 and 15<x<26. Performance improvements. Prototype 107 and 15<x<26. Performance improvements. Prototype
106 now supplies default value for relerr. 108 now supplies default value for relerr.
107 24 October 2012: Switch to continued-fraction expansion for 109 24 October 2012: Switch to continued-fraction expansion for
108 sufficiently large z, for performance reasons. 110 sufficiently large z, for performance reasons.
109 Also, avoid spurious overflow for |z| > 1e154. 111 Also, avoid spurious overflow for |z| > 1e154.
110 Set relerr argument to min(relerr,0.1). 112 Set relerr argument to min(relerr,0.1).
111 27 October 2012: Enhance accuracy in Re[w(z)] taken by itself, 113 27 October 2012: Enhance accuracy in Re[w(z)] taken by itself,
112 by switching to Alg. 916 in a region near 114 by switching to Alg. 916 in a region near
113 the real-z axis where continued fractions 115 the real-z axis where continued fractions
114 have poor relative accuracy in Re[w(z)]. Thanks 116 have poor relative accuracy in Re[w(z)]. Thanks
115 to M. Zaghloul for the tip. 117 to M. Zaghloul for the tip.
116 29 October 2012: Replace SLATEC-derived erfcx routine with 118 29 October 2012: Replace SLATEC-derived erfcx routine with
117 completely rewritten code by me, using a very 119 completely rewritten code by me, using a very
118 different algorithm which is much faster. 120 different algorithm which is much faster.
119 30 October 2012: Implemented special-case code for real z 121 30 October 2012: Implemented special-case code for real z
120 (where real part is exp(-x^2) and imag part is 122 (where real part is exp(-x^2) and imag part is
121 Dawson integral), using algorithm similar to erfx. 123 Dawson integral), using algorithm similar to erfx.
122 Export ImFaddeeva_w function to make Dawson's 124 Export ImFaddeeva_w function to make Dawson's
123 integral directly accessible. 125 integral directly accessible.
124 3 November 2012: Provide implementations of erf, erfc, erfcx, 126 3 November 2012: Provide implementations of erf, erfc, erfcx,
125 and Dawson functions in Faddeeva:: namespace, 127 and Dawson functions in Faddeeva:: namespace,
126 in addition to Faddeeva::w. Provide header 128 in addition to Faddeeva::w. Provide header
127 file Faddeeva.hh. 129 file Faddeeva.hh.
128 4 November 2012: Slightly faster erf for real arguments. 130 4 November 2012: Slightly faster erf for real arguments.
129 Updated MATLAB and Octave plugins. 131 Updated MATLAB and Octave plugins.
130 27 November 2012: Support compilation with either C++ or 132 27 November 2012: Support compilation with either C++ or
131 plain C (using C99 complex numbers). 133 plain C (using C99 complex numbers).
132 For real x, use standard-library erf(x) 134 For real x, use standard-library erf(x)
133 and erfc(x) if available (for C99 or C++11). 135 and erfc(x) if available (for C99 or C++11).
134 #include "config.h" if HAVE_CONFIG_H is #defined. 136 #include "config.h" if HAVE_CONFIG_H is #defined.
137 15 December 2012: Portability fixes (copysign, Inf/NaN creation),
138 use CMPLX/__builtin_complex if available in C,
139 slight accuracy improvements to erf and dawson
140 functions near the origin. Use gnulib functions
141 if GNULIB_NAMESPACE is defined.
135 */ 142 */
136 143
137 ///////////////////////////////////////////////////////////////////////// 144 /////////////////////////////////////////////////////////////////////////
138 /* If this file is compiled as a part of a larger project, 145 /* If this file is compiled as a part of a larger project,
139 support using an autoconf-style config.h header file 146 support using an autoconf-style config.h header file
143 #ifdef HAVE_CONFIG_H 150 #ifdef HAVE_CONFIG_H
144 # include "config.h" 151 # include "config.h"
145 #endif 152 #endif
146 153
147 ///////////////////////////////////////////////////////////////////////// 154 /////////////////////////////////////////////////////////////////////////
148 // basic math stuff, C++ & C99 compatibility 155 // macros to allow us to use either C++ or C (with C99 features)
149
150 #define Inf (1./0.) // infinity
151 #define NaN (0./0.) // NaN
152 156
153 #ifdef __cplusplus 157 #ifdef __cplusplus
154 158
155 # include "Faddeeva.hh" 159 # include "Faddeeva.hh"
156 160
157 # include <cfloat> 161 # include <cfloat>
158 # include <cmath> 162 # include <cmath>
163 # include <limits>
159 using namespace std; 164 using namespace std;
165
166 // use std::numeric_limits, since 1./0. and 0./0. fail with some compilers (MS)
167 # define Inf numeric_limits<double>::infinity()
168 # define NaN numeric_limits<double>::quiet_NaN()
160 169
161 typedef complex<double> cmplx; 170 typedef complex<double> cmplx;
162 171
163 // Use C-like complex syntax, since the C syntax is more restrictive 172 // Use C-like complex syntax, since the C syntax is more restrictive
164 # define cexp(z) exp(z) 173 # define cexp(z) exp(z)
165 # define creal(z) real(z) 174 # define creal(z) real(z)
166 # define cimag(z) imag(z) 175 # define cimag(z) imag(z)
167 # define cpolar(r,t) polar(r,t) 176 # define cpolar(r,t) polar(r,t)
168 177
169 # define CMPLX(a,b) cmplx(a,b) 178 # define C(a,b) cmplx(a,b)
170 179
171 # define FADDEEVA(name) Faddeeva::name 180 # define FADDEEVA(name) Faddeeva::name
172 # define FADDEEVA_RE(name) Faddeeva::name 181 # define FADDEEVA_RE(name) Faddeeva::name
173 182
174 // isnan/isinf were introduced in C++11 183 // isnan/isinf were introduced in C++11
175 # if (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF)) 184 # if (__cplusplus < 201103L) && (!defined(HAVE_ISNAN) || !defined(HAVE_ISINF))
176 static inline bool my_isnan(double x) { return x != x; } 185 static inline bool my_isnan(double x) { return x != x; }
177 # define isnan my_isnan 186 # define isnan my_isnan
178 static inline bool my_isinf(double x) { return 1/x == 0.; } 187 static inline bool my_isinf(double x) { return 1/x == 0.; }
179 # define isinf my_isinf 188 # define isinf my_isinf
189 # elif (__cplusplus >= 201103L)
190 // g++ gets confused between the C and C++ isnan/isinf functions
191 # define isnan std::isnan
192 # define isinf std::isinf
180 # endif 193 # endif
181 194
195 // copysign was introduced in C++11 (and is also in POSIX and C99)
196 # if defined(_WIN32) || defined(__WIN32__)
197 # define copysign _copysign // of course MS had to be different
198 # elif defined(GNULIB_NAMESPACE) // we are using using gnulib <cmath>
199 # define copysign GNULIB_NAMESPACE::copysign
200 # elif (__cplusplus < 201103L) && !defined(HAVE_COPYSIGN) && !defined(__linux__) && !(defined(__APPLE__) && defined(__MACH__)) && !defined(_AIX)
201 static inline double my_copysign(double x, double y) { return y<0 ? -x : x; }
202 # define copysign my_copysign
203 # endif
204
205 // If we are using the gnulib <cmath> (e.g. in the GNU Octave sources),
206 // gnulib generates a link warning if we use ::floor instead of gnulib::floor.
207 // This warning is completely innocuous because the only difference between
208 // gnulib::floor and the system ::floor (and only on ancient OSF systems)
209 // has to do with floor(-0), which doesn't occur in the usage below, but
210 // the Octave developers prefer that we silence the warning.
211 # ifdef GNULIB_NAMESPACE
212 # define floor GNULIB_NAMESPACE::floor
213 # endif
214
182 #else // !__cplusplus, i.e. pure C (requires C99 features) 215 #else // !__cplusplus, i.e. pure C (requires C99 features)
183 216
184 # include "Faddeeva.h" 217 # include "Faddeeva.h"
218
219 # define _GNU_SOURCE // enable GNU libc NAN extension if possible
185 220
186 # include <float.h> 221 # include <float.h>
187 # include <math.h> 222 # include <math.h>
188 223
189 typedef double complex cmplx; 224 typedef double complex cmplx;
190 225
191 # define CMPLX(a,b) ((a) + I*(b)) 226 # define FADDEEVA(name) Faddeeva_ ## name
227 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
228
229 /* Constructing complex numbers like 0+i*NaN is problematic in C99
230 without the C11 CMPLX macro, because 0.+I*NAN may give NaN+i*NAN if
231 I is a complex (rather than imaginary) constant. For some reason,
232 however, it works fine in (pre-4.7) gcc if I define Inf and NaN as
233 1/0 and 0/0 (and only if I compile with optimization -O1 or more),
234 but not if I use the INFINITY or NAN macros. */
235
236 /* __builtin_complex was introduced in gcc 4.7, but the C11 CMPLX macro
237 may not be defined unless we are using a recent (2012) version of
238 glibc and compile with -std=c11... note that icc lies about being
239 gcc and probably doesn't have this builtin(?), so exclude icc explicitly */
240 # if !defined(CMPLX) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !(defined(__ICC) || defined(__INTEL_COMPILER))
241 # define CMPLX(a,b) __builtin_complex((double) (a), (double) (b))
242 # endif
243
244 # ifdef CMPLX // C11
245 # define C(a,b) CMPLX(a,b)
246 # define Inf INFINITY // C99 infinity
247 # ifdef NAN // GNU libc extension
248 # define NaN NAN
249 # else
250 # define NaN (0./0.) // NaN
251 # endif
252 # else
253 # define C(a,b) ((a) + I*(b))
254 # define Inf (1./0.)
255 # define NaN (0./0.)
256 # endif
192 257
193 static inline cmplx cpolar(double r, double t) 258 static inline cmplx cpolar(double r, double t)
194 { 259 {
195 if (r == 0.0 && !isnan(t)) 260 if (r == 0.0 && !isnan(t))
196 return 0.0; 261 return 0.0;
197 else 262 else
198 return CMPLX(r * cos(t), r * sin(t)); 263 return C(r * cos(t), r * sin(t));
199 } 264 }
200 265
201 # define FADDEEVA(name) Faddeeva_ ## name 266 #endif // !__cplusplus, i.e. pure C (requires C99 features)
202 # define FADDEEVA_RE(name) Faddeeva_ ## name ## _re
203
204 #endif
205 267
206 ///////////////////////////////////////////////////////////////////////// 268 /////////////////////////////////////////////////////////////////////////
207 // Auxiliary routines to compute other special functions based on w(z) 269 // Auxiliary routines to compute other special functions based on w(z)
208 270
209 // compute erfcx(z) = exp(z^2) erfz(z) 271 // compute erfcx(z) = exp(z^2) erfz(z)
210 cmplx FADDEEVA(erfcx)(cmplx z, double relerr) 272 cmplx FADDEEVA(erfcx)(cmplx z, double relerr)
211 { 273 {
212 return FADDEEVA(w)(CMPLX(-cimag(z), creal(z)), relerr); 274 return FADDEEVA(w)(C(-cimag(z), creal(z)), relerr);
213 } 275 }
214 276
215 // compute the error function erf(x) 277 // compute the error function erf(x)
216 double FADDEEVA_RE(erf)(double x) 278 double FADDEEVA_RE(erf)(double x)
217 { 279 {
223 double mx2 = -x*x; 285 double mx2 = -x*x;
224 if (mx2 < -750) // underflow 286 if (mx2 < -750) // underflow
225 return (x >= 0 ? 1.0 : -1.0); 287 return (x >= 0 ? 1.0 : -1.0);
226 288
227 if (x >= 0) { 289 if (x >= 0) {
228 if (x < 5e-3) goto taylor; 290 if (x < 8e-2) goto taylor;
229 return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x); 291 return 1.0 - exp(mx2) * FADDEEVA_RE(erfcx)(x);
230 } 292 }
231 else { // x < 0 293 else { // x < 0
232 if (x > -5e-3) goto taylor; 294 if (x > -8e-2) goto taylor;
233 return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0; 295 return exp(mx2) * FADDEEVA_RE(erfcx)(-x) - 1.0;
234 } 296 }
235 297
236 // Use Taylor series for small |x|, to avoid cancellation inaccuracy 298 // Use Taylor series for small |x|, to avoid cancellation inaccuracy
237 // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - ...) 299 // erf(x) = 2/sqrt(pi) * x * (1 - x^2/3 + x^4/10 - x^6/42 + x^8/216 + ...)
238 taylor: 300 taylor:
239 return x * (1.1283791670955125739 301 return x * (1.1283791670955125739
240 + mx2 * (0.37612638903183752464 302 + mx2 * (0.37612638903183752464
241 + mx2 * 0.11283791670955125739)); 303 + mx2 * (0.11283791670955125739
304 + mx2 * (0.026866170645131251760
305 + mx2 * 0.0052239776254421878422))));
242 #endif 306 #endif
243 } 307 }
244 308
245 // compute the error function erf(z) 309 // compute the error function erf(z)
246 cmplx FADDEEVA(erf)(cmplx z, double relerr) 310 cmplx FADDEEVA(erf)(cmplx z, double relerr)
247 { 311 {
248 double x = creal(z), y = cimag(z); 312 double x = creal(z), y = cimag(z);
249 313
250 if (y == 0) 314 if (y == 0)
251 return CMPLX(FADDEEVA_RE(erf)(x), 315 return C(FADDEEVA_RE(erf)(x),
252 y); // preserve sign of 0 316 y); // preserve sign of 0
253 if (x == 0) // handle separately for speed & handling of y = Inf or NaN 317 if (x == 0) // handle separately for speed & handling of y = Inf or NaN
254 return CMPLX(x, // preserve sign of 0 318 return C(x, // preserve sign of 0
255 /* handle y -> Inf limit manually, since 319 /* handle y -> Inf limit manually, since
256 exp(y^2) -> Inf but Im[w(y)] -> 0, so 320 exp(y^2) -> Inf but Im[w(y)] -> 0, so
257 IEEE will give us a NaN when it should be Inf */ 321 IEEE will give us a NaN when it should be Inf */
258 y*y > 720 ? (y > 0 ? Inf : -Inf) 322 y*y > 720 ? (y > 0 ? Inf : -Inf)
259 : exp(y*y) * FADDEEVA(w_im)(y)); 323 : exp(y*y) * FADDEEVA(w_im)(y));
260 324
261 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 325 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
262 double mIm_z2 = -2*x*y; // Im(-z^2) 326 double mIm_z2 = -2*x*y; // Im(-z^2)
263 if (mRe_z2 < -750) // underflow 327 if (mRe_z2 < -750) // underflow
264 return (x >= 0 ? 1.0 : -1.0); 328 return (x >= 0 ? 1.0 : -1.0);
265 329
266 /* Handle positive and negative x via different formulas, 330 /* Handle positive and negative x via different formulas,
267 using the mirror symmetries of w, to avoid overflow/underflow 331 using the mirror symmetries of w, to avoid overflow/underflow
268 problems from multiplying exponentially large and small quantities. */ 332 problems from multiplying exponentially large and small quantities. */
269 if (x >= 0) { 333 if (x >= 0) {
270 if (x < 5e-3) { 334 if (x < 8e-2) {
271 if (fabs(y) < 5e-3) 335 if (fabs(y) < 1e-2)
272 goto taylor; 336 goto taylor;
273 else if (fabs(mIm_z2) < 5e-3) 337 else if (fabs(mIm_z2) < 5e-3 && x < 5e-3)
274 goto taylor_erfi; 338 goto taylor_erfi;
275 } 339 }
276 /* don't use complex exp function, since that will produce spurious NaN 340 /* don't use complex exp function, since that will produce spurious NaN
277 values when multiplying w in an overflow situation. */ 341 values when multiplying w in an overflow situation. */
278 return 1.0 - exp(mRe_z2) * 342 return 1.0 - exp(mRe_z2) *
279 (CMPLX(cos(mIm_z2), sin(mIm_z2)) 343 (C(cos(mIm_z2), sin(mIm_z2))
280 * FADDEEVA(w)(CMPLX(-y,x), relerr)); 344 * FADDEEVA(w)(C(-y,x), relerr));
281 } 345 }
282 else { // x < 0 346 else { // x < 0
283 if (x > -5e-3) { // duplicate from above to avoid fabs(x) call 347 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call
284 if (fabs(y) < 5e-3) 348 if (fabs(y) < 1e-2)
285 goto taylor; 349 goto taylor;
286 else if (fabs(mIm_z2) < 5e-3) 350 else if (fabs(mIm_z2) < 5e-3 && x > -5e-3)
287 goto taylor_erfi; 351 goto taylor_erfi;
288 } 352 }
289 else if (isnan(x)) 353 else if (isnan(x))
290 return CMPLX(NaN, y == 0 ? 0 : NaN); 354 return C(NaN, y == 0 ? 0 : NaN);
291 /* don't use complex exp function, since that will produce spurious NaN 355 /* don't use complex exp function, since that will produce spurious NaN
292 values when multiplying w in an overflow situation. */ 356 values when multiplying w in an overflow situation. */
293 return exp(mRe_z2) * 357 return exp(mRe_z2) *
294 (CMPLX(cos(mIm_z2), sin(mIm_z2)) 358 (C(cos(mIm_z2), sin(mIm_z2))
295 * FADDEEVA(w)(CMPLX(y,-x), relerr)) - 1.0; 359 * FADDEEVA(w)(C(y,-x), relerr)) - 1.0;
296 } 360 }
297 361
298 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 362 // Use Taylor series for small |z|, to avoid cancellation inaccuracy
299 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - ...) 363 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...)
300 taylor: 364 taylor:
301 { 365 {
302 cmplx mz2 = CMPLX(mRe_z2, mIm_z2); // -z^2 366 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
303 return z * (1.1283791670955125739 367 return z * (1.1283791670955125739
304 + mz2 * (0.37612638903183752464 368 + mz2 * (0.37612638903183752464
305 + mz2 * 0.11283791670955125739)); 369 + mz2 * (0.11283791670955125739
370 + mz2 * (0.026866170645131251760
371 + mz2 * 0.0052239776254421878422))));
306 } 372 }
307 373
308 /* for small |x| and small |xy|, 374 /* for small |x| and small |xy|,
309 use Taylor series to avoid cancellation inaccuracy: 375 use Taylor series to avoid cancellation inaccuracy:
310 erf(x+iy) = erf(iy) 376 erf(x+iy) = erf(iy)
311 + 2*exp(y^2)/sqrt(pi) * 377 + 2*exp(y^2)/sqrt(pi) *
312 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 378 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ...
313 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] 379 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ]
314 where: 380 where:
315 erf(iy) = exp(y^2) * Im[w(y)] 381 erf(iy) = exp(y^2) * Im[w(y)]
316 */ 382 */
317 taylor_erfi: 383 taylor_erfi:
318 { 384 {
319 double x2 = x*x, y2 = y*y; 385 double x2 = x*x, y2 = y*y;
320 double expy2 = exp(y2); 386 double expy2 = exp(y2);
321 return CMPLX 387 return C
322 (expy2 * x * (1.1283791670955125739 388 (expy2 * x * (1.1283791670955125739
323 - x2 * (0.37612638903183752464 389 - x2 * (0.37612638903183752464
324 + 0.75225277806367504925*y2) 390 + 0.75225277806367504925*y2)
325 + x2*x2 * (0.11283791670955125739 391 + x2*x2 * (0.11283791670955125739
326 + y2 * (0.45135166683820502956 392 + y2 * (0.45135166683820502956
327 + 0.15045055561273500986*y2))), 393 + 0.15045055561273500986*y2))),
328 expy2 * (FADDEEVA(w_im)(y) 394 expy2 * (FADDEEVA(w_im)(y)
329 - x2*y * (1.1283791670955125739 395 - x2*y * (1.1283791670955125739
330 - x2 * (0.56418958354775628695 396 - x2 * (0.56418958354775628695
331 + 0.37612638903183752464*y2)))); 397 + 0.37612638903183752464*y2))));
332 } 398 }
333 } 399 }
334 400
335 // erfi(z) = -i erf(iz) 401 // erfi(z) = -i erf(iz)
336 cmplx FADDEEVA(erfi)(cmplx z, double relerr) 402 cmplx FADDEEVA(erfi)(cmplx z, double relerr)
337 { 403 {
338 cmplx e = FADDEEVA(erf)(CMPLX(-cimag(z),creal(z)), relerr); 404 cmplx e = FADDEEVA(erf)(C(-cimag(z),creal(z)), relerr);
339 return CMPLX(cimag(e), -creal(e)); 405 return C(cimag(e), -creal(e));
340 } 406 }
341 407
342 // erfi(x) = -i erf(ix) 408 // erfi(x) = -i erf(ix)
343 double FADDEEVA_RE(erfi)(double x) 409 double FADDEEVA_RE(erfi)(double x)
344 { 410 {
365 cmplx FADDEEVA(erfc)(cmplx z, double relerr) 431 cmplx FADDEEVA(erfc)(cmplx z, double relerr)
366 { 432 {
367 double x = creal(z), y = cimag(z); 433 double x = creal(z), y = cimag(z);
368 434
369 if (x == 0.) 435 if (x == 0.)
370 return CMPLX(1, 436 return C(1,
371 /* handle y -> Inf limit manually, since 437 /* handle y -> Inf limit manually, since
372 exp(y^2) -> Inf but Im[w(y)] -> 0, so 438 exp(y^2) -> Inf but Im[w(y)] -> 0, so
373 IEEE will give us a NaN when it should be Inf */ 439 IEEE will give us a NaN when it should be Inf */
374 y*y > 720 ? (y > 0 ? -Inf : Inf) 440 y*y > 720 ? (y > 0 ? -Inf : Inf)
375 : -exp(y*y) * FADDEEVA(w_im)(y)); 441 : -exp(y*y) * FADDEEVA(w_im)(y));
376 if (y == 0.) { 442 if (y == 0.) {
377 if (x*x > 750) // underflow 443 if (x*x > 750) // underflow
378 return CMPLX(x >= 0 ? 0.0 : 2.0, 444 return C(x >= 0 ? 0.0 : 2.0,
379 -y); // preserve sign of 0 445 -y); // preserve sign of 0
380 return CMPLX(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x) 446 return C(x >= 0 ? exp(-x*x) * FADDEEVA_RE(erfcx)(x)
381 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x), 447 : 2. - exp(-x*x) * FADDEEVA_RE(erfcx)(-x),
382 -y); // preserve sign of zero 448 -y); // preserve sign of zero
383 } 449 }
384 450
385 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 451 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
386 double mIm_z2 = -2*x*y; // Im(-z^2) 452 double mIm_z2 = -2*x*y; // Im(-z^2)
387 if (mRe_z2 < -750) // underflow 453 if (mRe_z2 < -750) // underflow
388 return (x >= 0 ? 0.0 : 2.0); 454 return (x >= 0 ? 0.0 : 2.0);
389 455
390 if (x >= 0) 456 if (x >= 0)
391 return cexp(CMPLX(mRe_z2, mIm_z2)) 457 return cexp(C(mRe_z2, mIm_z2))
392 * FADDEEVA(w)(CMPLX(-y,x), relerr); 458 * FADDEEVA(w)(C(-y,x), relerr);
393 else 459 else
394 return 2.0 - cexp(CMPLX(mRe_z2, mIm_z2)) 460 return 2.0 - cexp(C(mRe_z2, mIm_z2))
395 * FADDEEVA(w)(CMPLX(y,-x), relerr); 461 * FADDEEVA(w)(C(y,-x), relerr);
396 } 462 }
397 463
398 // compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x) 464 // compute Dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x)
399 double FADDEEVA_RE(Dawson)(double x) 465 double FADDEEVA_RE(Dawson)(double x)
400 { 466 {
408 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 474 const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2
409 double x = creal(z), y = cimag(z); 475 double x = creal(z), y = cimag(z);
410 476
411 // handle axes separately for speed & proper handling of x or y = Inf or NaN 477 // handle axes separately for speed & proper handling of x or y = Inf or NaN
412 if (y == 0) 478 if (y == 0)
413 return CMPLX(spi2 * FADDEEVA(w_im)(x), 479 return C(spi2 * FADDEEVA(w_im)(x),
414 -y); // preserve sign of 0 480 -y); // preserve sign of 0
415 if (x == 0) { 481 if (x == 0) {
416 double y2 = y*y; 482 double y2 = y*y;
417 if (y2 < 2.5e-5) { // Taylor expansion 483 if (y2 < 2.5e-5) { // Taylor expansion
418 return CMPLX(x, // preserve sign of 0 484 return C(x, // preserve sign of 0
419 y * (1. 485 y * (1.
420 + y2 * (0.6666666666666666666666666666666666666667 486 + y2 * (0.6666666666666666666666666666666666666667
421 + y2 * 0.2666666666666666666666666666666666666667))); 487 + y2 * 0.26666666666666666666666666666666666667)));
422 } 488 }
423 return CMPLX(x, // preserve sign of 0 489 return C(x, // preserve sign of 0
424 spi2 * (y >= 0 490 spi2 * (y >= 0
425 ? exp(y2) - FADDEEVA_RE(erfcx)(y) 491 ? exp(y2) - FADDEEVA_RE(erfcx)(y)
426 : FADDEEVA_RE(erfcx)(-y) - exp(y2))); 492 : FADDEEVA_RE(erfcx)(-y) - exp(y2)));
427 } 493 }
428 494
429 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 495 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow
430 double mIm_z2 = -2*x*y; // Im(-z^2) 496 double mIm_z2 = -2*x*y; // Im(-z^2)
431 cmplx mz2 = CMPLX(mRe_z2, mIm_z2); // -z^2 497 cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2
432 498
433 /* Handle positive and negative x via different formulas, 499 /* Handle positive and negative x via different formulas,
434 using the mirror symmetries of w, to avoid overflow/underflow 500 using the mirror symmetries of w, to avoid overflow/underflow
435 problems from multiplying exponentially large and small quantities. */ 501 problems from multiplying exponentially large and small quantities. */
436 if (y >= 0) { 502 if (y >= 0) {
437 if (y < 5e-3) { 503 if (y < 5e-3) {
438 if (fabs(x) < 5e-3) 504 if (fabs(x) < 5e-3)
439 goto taylor; 505 goto taylor;
440 else if (fabs(mIm_z2) < 5e-3) 506 else if (fabs(mIm_z2) < 5e-3)
441 goto taylor_realaxis; 507 goto taylor_realaxis;
442 } 508 }
443 cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr); 509 cmplx res = cexp(mz2) - FADDEEVA(w)(z, relerr);
444 return spi2 * CMPLX(-cimag(res), creal(res)); 510 return spi2 * C(-cimag(res), creal(res));
445 } 511 }
446 else { // y < 0 512 else { // y < 0
447 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call 513 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call
448 if (fabs(x) < 5e-3) 514 if (fabs(x) < 5e-3)
449 goto taylor; 515 goto taylor;
450 else if (fabs(mIm_z2) < 5e-3) 516 else if (fabs(mIm_z2) < 5e-3)
451 goto taylor_realaxis; 517 goto taylor_realaxis;
452 } 518 }
453 else if (isnan(y)) 519 else if (isnan(y))
454 return CMPLX(x == 0 ? 0 : NaN, NaN); 520 return C(x == 0 ? 0 : NaN, NaN);
455 cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2); 521 cmplx res = FADDEEVA(w)(-z, relerr) - cexp(mz2);
456 return spi2 * CMPLX(-cimag(res), creal(res)); 522 return spi2 * C(-cimag(res), creal(res));
457 } 523 }
458 524
459 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 525 // Use Taylor series for small |z|, to avoid cancellation inaccuracy
460 // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... 526 // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ...
461 taylor: 527 taylor:
462 return z * (1. 528 return z * (1.
463 + mz2 * (0.6666666666666666666666666666666666666667 529 + mz2 * (0.6666666666666666666666666666666666666667
464 + mz2 * 0.2666666666666666666666666666666666666667)); 530 + mz2 * 0.2666666666666666666666666666666666666667));
465 531
466 /* for small |y| and small |xy|, 532 /* for small |y| and small |xy|,
467 use Taylor series to avoid cancellation inaccuracy: 533 use Taylor series to avoid cancellation inaccuracy:
468 dawson(x + iy) 534 dawson(x + iy)
469 = D + y^2 (D + x - 2Dx^2) 535 = D + y^2 (D + x - 2Dx^2)
500 { 566 {
501 double x2 = x*x; 567 double x2 = x*x;
502 if (x2 > 1600) { // |x| > 40 568 if (x2 > 1600) { // |x| > 40
503 double y2 = y*y; 569 double y2 = y*y;
504 if (x2 > 25e14) {// |x| > 5e7 570 if (x2 > 25e14) {// |x| > 5e7
505 double xy2 = (x*y)*(x*y); 571 double xy2 = (x*y)*(x*y);
506 return CMPLX((0.5 + y2 * (0.5 + 0.25*y2 572 return C((0.5 + y2 * (0.5 + 0.25*y2
507 - 0.16666666666666666667*xy2)) / x, 573 - 0.16666666666666666667*xy2)) / x,
508 y * (-1 + y2 * (-0.66666666666666666667 574 y * (-1 + y2 * (-0.66666666666666666667
509 + 0.13333333333333333333*xy2 575 + 0.13333333333333333333*xy2
510 - 0.26666666666666666667*y2)) 576 - 0.26666666666666666667*y2))
511 / (2*x2 - 1)); 577 / (2*x2 - 1));
512 } 578 }
513 return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) * 579 return (1. / (-15 + x2*(90 + x2*(-60 + 8*x2)))) *
514 CMPLX(x * (33 + x2 * (-28 + 4*x2) 580 C(x * (33 + x2 * (-28 + 4*x2)
515 + y2 * (18 - 4*x2 + 4*y2)), 581 + y2 * (18 - 4*x2 + 4*y2)),
516 y * (-15 + x2 * (24 - 4*x2) 582 y * (-15 + x2 * (24 - 4*x2)
517 + y2 * (4*x2 - 10 - 4*y2))); 583 + y2 * (4*x2 - 10 - 4*y2)));
518 } 584 }
519 else { 585 else {
520 double D = spi2 * FADDEEVA(w_im)(x); 586 double D = spi2 * FADDEEVA(w_im)(x);
521 double x2 = x*x, y2 = y*y; 587 double x2 = x*x, y2 = y*y;
522 return CMPLX 588 return C
523 (D + y2 * (D + x - 2*D*x2) 589 (D + y2 * (D + x - 2*D*x2)
524 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2)) 590 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2))
525 + x * (0.83333333333333333333 591 + x * (0.83333333333333333333
526 - 0.33333333333333333333 * x2)), 592 - 0.33333333333333333333 * x2)),
527 y * (1 - 2*D*x 593 y * (1 - 2*D*x
528 + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2)) 594 + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2))
529 + y2*y2 * (0.26666666666666666667 - 595 + y2*y2 * (0.26666666666666666667 -
530 x2 * (0.6 - 0.13333333333333333333 * x2) 596 x2 * (0.6 - 0.13333333333333333333 * x2)
531 - D*x * (1 - x2 * (1.3333333333333333333 597 - D*x * (1 - x2 * (1.3333333333333333333
532 - 0.26666666666666666667 * x2))))); 598 - 0.26666666666666666667 * x2)))));
533 } 599 }
534 } 600 }
535 } 601 }
536 602
537 ///////////////////////////////////////////////////////////////////////// 603 /////////////////////////////////////////////////////////////////////////
543 } 609 }
544 610
545 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2 611 // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2
546 static inline double sinh_taylor(double x) { 612 static inline double sinh_taylor(double x) {
547 return x * (1 + (x*x) * (0.1666666666666666666667 613 return x * (1 + (x*x) * (0.1666666666666666666667
548 + 0.00833333333333333333333 * (x*x))); 614 + 0.00833333333333333333333 * (x*x)));
549 } 615 }
550 616
551 static inline double sqr(double x) { return x*x; } 617 static inline double sqr(double x) { return x*x; }
552 618
553 // precomputed table of expa2n2[n-1] = exp(-a2*n*n) 619 // precomputed table of expa2n2[n-1] = exp(-a2*n*n)
610 ///////////////////////////////////////////////////////////////////////// 676 /////////////////////////////////////////////////////////////////////////
611 677
612 cmplx FADDEEVA(w)(cmplx z, double relerr) 678 cmplx FADDEEVA(w)(cmplx z, double relerr)
613 { 679 {
614 if (creal(z) == 0.0) 680 if (creal(z) == 0.0)
615 return CMPLX(FADDEEVA_RE(erfcx)(cimag(z)), 681 return C(FADDEEVA_RE(erfcx)(cimag(z)),
616 creal(z)); // give correct sign of 0 in cimag(w) 682 creal(z)); // give correct sign of 0 in cimag(w)
617 else if (cimag(z) == 0) 683 else if (cimag(z) == 0)
618 return CMPLX(exp(-sqr(creal(z))), 684 return C(exp(-sqr(creal(z))),
619 FADDEEVA(w_im)(creal(z))); 685 FADDEEVA(w_im)(creal(z)));
620 686
621 double a, a2, c; 687 double a, a2, c;
622 if (relerr <= DBL_EPSILON) { 688 if (relerr <= DBL_EPSILON) {
623 relerr = DBL_EPSILON; 689 relerr = DBL_EPSILON;
624 a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5)) 690 a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5))
641 707
642 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z| 708 #define USE_CONTINUED_FRACTION 1 // 1 to use continued fraction for large |z|
643 709
644 #if USE_CONTINUED_FRACTION 710 #if USE_CONTINUED_FRACTION
645 if (ya > 7 || (x > 6 // continued fraction is faster 711 if (ya > 7 || (x > 6 // continued fraction is faster
646 /* As pointed out by M. Zaghloul, the continued 712 /* As pointed out by M. Zaghloul, the continued
647 fraction seems to give a large relative error in 713 fraction seems to give a large relative error in
648 Re w(z) for |x| ~ 6 and small |y|, so use 714 Re w(z) for |x| ~ 6 and small |y|, so use
649 algorithm 816 in this region: */ 715 algorithm 816 in this region: */
650 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) { 716 && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) {
651 717
652 /* Poppe & Wijers suggest using a number of terms 718 /* Poppe & Wijers suggest using a number of terms
653 nu = 3 + 1442 / (26*rho + 77) 719 nu = 3 + 1442 / (26*rho + 77)
654 where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4. 720 where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4.
655 (They only use this expansion for rho >= 1, but rho a little less 721 (They only use this expansion for rho >= 1, but rho a little less
661 I also separate the regions where nu == 2 and nu == 1. */ 727 I also separate the regions where nu == 2 and nu == 1. */
662 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 728 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
663 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 729 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
664 if (x + ya > 4000) { // nu <= 2 730 if (x + ya > 4000) { // nu <= 2
665 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z 731 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z
666 // scale to avoid overflow 732 // scale to avoid overflow
667 if (x > ya) { 733 if (x > ya) {
668 double yax = ya / xs; 734 double yax = ya / xs;
669 double denom = ispi / (xs + yax*ya); 735 double denom = ispi / (xs + yax*ya);
670 ret = CMPLX(denom*yax, denom); 736 ret = C(denom*yax, denom);
671 } 737 }
672 else if (isinf(ya)) 738 else if (isinf(ya))
673 return ((isnan(x) || y < 0) 739 return ((isnan(x) || y < 0)
674 ? CMPLX(NaN,NaN) : CMPLX(0,0)); 740 ? C(NaN,NaN) : C(0,0));
675 else { 741 else {
676 double xya = xs / ya; 742 double xya = xs / ya;
677 double denom = ispi / (xya*xs + ya); 743 double denom = ispi / (xya*xs + ya);
678 ret = CMPLX(denom, denom*xya); 744 ret = C(denom, denom*xya);
679 } 745 }
680 } 746 }
681 else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5) 747 else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5)
682 double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya; 748 double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya;
683 double denom = ispi / (dr*dr + di*di); 749 double denom = ispi / (dr*dr + di*di);
684 ret = CMPLX(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di)); 750 ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di));
685 } 751 }
686 } 752 }
687 else { // compute nu(z) estimate and do general continued fraction 753 else { // compute nu(z) estimate and do general continued fraction
688 const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit 754 const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit
689 double nu = floor(c0 + c1 / (c2*x + c3*ya + c4)); 755 double nu = floor(c0 + c1 / (c2*x + c3*ya + c4));
690 double wr = xs, wi = ya; 756 double wr = xs, wi = ya;
691 for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) { 757 for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) {
692 // w <- z - nu/w: 758 // w <- z - nu/w:
693 double denom = nu / (wr*wr + wi*wi); 759 double denom = nu / (wr*wr + wi*wi);
694 wr = xs - wr * denom; 760 wr = xs - wr * denom;
695 wi = ya + wi * denom; 761 wi = ya + wi * denom;
696 } 762 }
697 { // w(z) = i/sqrt(pi) / w: 763 { // w(z) = i/sqrt(pi) / w:
698 double denom = ispi / (wr*wr + wi*wi); 764 double denom = ispi / (wr*wr + wi*wi);
699 ret = CMPLX(denom*wi, denom*wr); 765 ret = C(denom*wi, denom*wr);
700 } 766 }
701 } 767 }
702 if (y < 0) { 768 if (y < 0) {
703 // use w(z) = 2.0*exp(-z*z) - w(-z), 769 // use w(z) = 2.0*exp(-z*z) - w(-z),
704 // but be careful of overflow in exp(-z*z) 770 // but be careful of overflow in exp(-z*z)
705 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) 771 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
706 return 2.0*cexp(CMPLX((ya-xs)*(xs+ya), 2*xs*y)) - ret; 772 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
707 } 773 }
708 else 774 else
709 return ret; 775 return ret;
710 } 776 }
711 #else // !USE_CONTINUED_FRACTION 777 #else // !USE_CONTINUED_FRACTION
714 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 780 double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0
715 // scale to avoid overflow 781 // scale to avoid overflow
716 if (x > ya) { 782 if (x > ya) {
717 double yax = ya / xs; 783 double yax = ya / xs;
718 double denom = ispi / (xs + yax*ya); 784 double denom = ispi / (xs + yax*ya);
719 ret = CMPLX(denom*yax, denom); 785 ret = C(denom*yax, denom);
720 } 786 }
721 else { 787 else {
722 double xya = xs / ya; 788 double xya = xs / ya;
723 double denom = ispi / (xya*xs + ya); 789 double denom = ispi / (xya*xs + ya);
724 ret = CMPLX(denom, denom*xya); 790 ret = C(denom, denom*xya);
725 } 791 }
726 if (y < 0) { 792 if (y < 0) {
727 // use w(z) = 2.0*exp(-z*z) - w(-z), 793 // use w(z) = 2.0*exp(-z*z) - w(-z),
728 // but be careful of overflow in exp(-z*z) 794 // but be careful of overflow in exp(-z*z)
729 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) 795 // = exp(-(xs*xs-ya*ya) -2*i*xs*ya)
730 return 2.0*cexp(CMPLX((ya-xs)*(xs+ya), 2*xs*y)) - ret; 796 return 2.0*cexp(C((ya-xs)*(xs+ya), 2*xs*y)) - ret;
731 } 797 }
732 else 798 else
733 return ret; 799 return ret;
734 } 800 }
735 #endif // !USE_CONTINUED_FRACTION 801 #endif // !USE_CONTINUED_FRACTION
749 else if (x < 10) { 815 else if (x < 10) {
750 double prod2ax = 1, prodm2ax = 1; 816 double prod2ax = 1, prodm2ax = 1;
751 double expx2; 817 double expx2;
752 818
753 if (isnan(y)) 819 if (isnan(y))
754 return CMPLX(y,y); 820 return C(y,y);
755 821
756 /* Somewhat ugly copy-and-paste duplication here, but I see significant 822 /* Somewhat ugly copy-and-paste duplication here, but I see significant
757 speedups from using the special-case code with the precomputed 823 speedups from using the special-case code with the precomputed
758 exponential, and the x < 5e-4 special case is needed for accuracy. */ 824 exponential, and the x < 5e-4 special case is needed for accuracy. */
759 825
760 if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table 826 if (relerr == DBL_EPSILON) { // use precomputed exp(-a2*(n*n)) table
761 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 827 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
762 const double x2 = x*x; 828 const double x2 = x*x;
763 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor 829 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
764 // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision 830 // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision
765 const double ax2 = 1.036642960860171859744*x; // 2*a*x 831 const double ax2 = 1.036642960860171859744*x; // 2*a*x
766 const double exp2ax = 832 const double exp2ax =
767 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2)); 833 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2));
768 const double expm2ax = 834 const double expm2ax =
769 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2)); 835 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2));
770 for (int n = 1; 1; ++n) { 836 for (int n = 1; 1; ++n) {
771 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); 837 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
772 prod2ax *= exp2ax; 838 prod2ax *= exp2ax;
773 prodm2ax *= expm2ax; 839 prodm2ax *= expm2ax;
774 sum1 += coef; 840 sum1 += coef;
775 sum2 += coef * prodm2ax; 841 sum2 += coef * prodm2ax;
776 sum3 += coef * prod2ax; 842 sum3 += coef * prod2ax;
777 843
778 // really = sum5 - sum4 844 // really = sum5 - sum4
779 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x); 845 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
780 846
781 // test convergence via sum3 847 // test convergence via sum3
782 if (coef * prod2ax < relerr * sum3) break; 848 if (coef * prod2ax < relerr * sum3) break;
783 } 849 }
784 } 850 }
785 else { // x > 5e-4, compute sum4 and sum5 separately 851 else { // x > 5e-4, compute sum4 and sum5 separately
786 expx2 = exp(-x*x); 852 expx2 = exp(-x*x);
787 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax; 853 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
788 for (int n = 1; 1; ++n) { 854 for (int n = 1; 1; ++n) {
789 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); 855 const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y);
790 prod2ax *= exp2ax; 856 prod2ax *= exp2ax;
791 prodm2ax *= expm2ax; 857 prodm2ax *= expm2ax;
792 sum1 += coef; 858 sum1 += coef;
793 sum2 += coef * prodm2ax; 859 sum2 += coef * prodm2ax;
794 sum4 += (coef * prodm2ax) * (a*n); 860 sum4 += (coef * prodm2ax) * (a*n);
795 sum3 += coef * prod2ax; 861 sum3 += coef * prod2ax;
796 sum5 += (coef * prod2ax) * (a*n); 862 sum5 += (coef * prod2ax) * (a*n);
797 // test convergence via sum5, since this sum has the slowest decay 863 // test convergence via sum5, since this sum has the slowest decay
798 if ((coef * prod2ax) * (a*n) < relerr * sum5) break; 864 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
799 } 865 }
800 } 866 }
801 } 867 }
802 else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly 868 else { // relerr != DBL_EPSILON, compute exp(-a2*(n*n)) on the fly
803 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax; 869 const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax;
804 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 870 if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4
805 const double x2 = x*x; 871 const double x2 = x*x;
806 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor 872 expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor
807 for (int n = 1; 1; ++n) { 873 for (int n = 1; 1; ++n) {
808 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y); 874 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
809 prod2ax *= exp2ax; 875 prod2ax *= exp2ax;
810 prodm2ax *= expm2ax; 876 prodm2ax *= expm2ax;
811 sum1 += coef; 877 sum1 += coef;
812 sum2 += coef * prodm2ax; 878 sum2 += coef * prodm2ax;
813 sum3 += coef * prod2ax; 879 sum3 += coef * prod2ax;
814 880
815 // really = sum5 - sum4 881 // really = sum5 - sum4
816 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x); 882 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x);
817 883
818 // test convergence via sum3 884 // test convergence via sum3
819 if (coef * prod2ax < relerr * sum3) break; 885 if (coef * prod2ax < relerr * sum3) break;
820 } 886 }
821 } 887 }
822 else { // x > 5e-4, compute sum4 and sum5 separately 888 else { // x > 5e-4, compute sum4 and sum5 separately
823 expx2 = exp(-x*x); 889 expx2 = exp(-x*x);
824 for (int n = 1; 1; ++n) { 890 for (int n = 1; 1; ++n) {
825 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y); 891 const double coef = exp(-a2*(n*n)) * expx2 / (a2*(n*n) + y*y);
826 prod2ax *= exp2ax; 892 prod2ax *= exp2ax;
827 prodm2ax *= expm2ax; 893 prodm2ax *= expm2ax;
828 sum1 += coef; 894 sum1 += coef;
829 sum2 += coef * prodm2ax; 895 sum2 += coef * prodm2ax;
830 sum4 += (coef * prodm2ax) * (a*n); 896 sum4 += (coef * prodm2ax) * (a*n);
831 sum3 += coef * prod2ax; 897 sum3 += coef * prod2ax;
832 sum5 += (coef * prod2ax) * (a*n); 898 sum5 += (coef * prod2ax) * (a*n);
833 // test convergence via sum5, since this sum has the slowest decay 899 // test convergence via sum5, since this sum has the slowest decay
834 if ((coef * prod2ax) * (a*n) < relerr * sum5) break; 900 if ((coef * prod2ax) * (a*n) < relerr * sum5) break;
835 } 901 }
836 } 902 }
837 } 903 }
838 const double expx2erfcxy = // avoid spurious overflow for large negative y 904 const double expx2erfcxy = // avoid spurious overflow for large negative y
839 y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision 905 y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision
840 ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x); 906 ? expx2*FADDEEVA_RE(erfcx)(y) : 2*exp(y*y-x*x);
841 if (y > 5) { // imaginary terms cancel 907 if (y > 5) { // imaginary terms cancel
842 const double sinxy = sin(x*y); 908 const double sinxy = sin(x*y);
843 ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y) 909 ret = (expx2erfcxy - c*y*sum1) * cos(2*x*y)
844 + (c*x*expx2) * sinxy * sinc(x*y, sinxy); 910 + (c*x*expx2) * sinxy * sinc(x*y, sinxy);
845 } 911 }
846 else { 912 else {
847 double xs = creal(z); 913 double xs = creal(z);
848 const double sinxy = sin(xs*y); 914 const double sinxy = sin(xs*y);
849 const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y); 915 const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y);
850 const double coef1 = expx2erfcxy - c*y*sum1; 916 const double coef1 = expx2erfcxy - c*y*sum1;
851 const double coef2 = c*xs*expx2; 917 const double coef2 = c*xs*expx2;
852 ret = CMPLX(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy), 918 ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy),
853 coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy); 919 coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy);
854 } 920 }
855 } 921 }
856 else { // x large: only sum3 & sum5 contribute (see above note) 922 else { // x large: only sum3 & sum5 contribute (see above note)
857 if (isnan(x)) 923 if (isnan(x))
858 return CMPLX(x,x); 924 return C(x,x);
859 if (isnan(y)) 925 if (isnan(y))
860 return CMPLX(y,y); 926 return C(y,y);
861 927
862 #if USE_CONTINUED_FRACTION 928 #if USE_CONTINUED_FRACTION
863 ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term 929 ret = exp(-x*x); // |y| < 1e-10, so we only need exp(-x*x) term
864 #else 930 #else
865 if (y < 0) { 931 if (y < 0) {
866 /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so 932 /* erfcx(y) ~ 2*exp(y*y) + (< 1) if y < 0, so
867 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible 933 erfcx(y)*exp(-x*x) ~ 2*exp(y*y-x*x) term may not be negligible
868 if y*y - x*x > -36 or so. So, compute this term just in case. 934 if y*y - x*x > -36 or so. So, compute this term just in case.
869 We also need the -exp(-x*x) term to compute Re[w] accurately 935 We also need the -exp(-x*x) term to compute Re[w] accurately
870 in the case where y is very small. */ 936 in the case where y is very small. */
871 ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y); 937 ret = cpolar(2*exp(y*y-x*x) - exp(-x*x), -2*creal(z)*y);
872 } 938 }
873 else 939 else
874 ret = exp(-x*x); // not negligible in real part if y very small 940 ret = exp(-x*x); // not negligible in real part if y very small
875 #endif 941 #endif
897 sum5 += a * np * tp; 963 sum5 += a * np * tp;
898 if (a * np * tp < relerr * sum5) goto finish; 964 if (a * np * tp < relerr * sum5) goto finish;
899 } 965 }
900 } 966 }
901 finish: 967 finish:
902 return ret + CMPLX((0.5*c)*y*(sum2+sum3), 968 return ret + C((0.5*c)*y*(sum2+sum3),
903 (0.5*c)*copysign(sum5-sum4, creal(z))); 969 (0.5*c)*copysign(sum5-sum4, creal(z)));
904 } 970 }
905 971
906 ///////////////////////////////////////////////////////////////////////// 972 /////////////////////////////////////////////////////////////////////////
907 973
908 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by 974 /* erfcx(x) = exp(x^2) erfc(x) function, for real x, written by
916 Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations, 982 Second, for 0 <= x <= 50, it uses Chebyshev polynomial approximations,
917 but with two twists: 983 but with two twists:
918 984
919 a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation, 985 a) It maps x to y = 4 / (4+x) in [0,1]. This simple transformation,
920 inspired by a similar transformation in the octave-forge/specfun 986 inspired by a similar transformation in the octave-forge/specfun
921 erfcx by Soren Hauberg, results in much faster Chebyshev convergence 987 erfcx by Soren Hauberg, results in much faster Chebyshev convergence
922 than other simple transformations I have examined. 988 than other simple transformations I have examined.
923 989
924 b) Instead of using a single Chebyshev polynomial for the entire 990 b) Instead of using a single Chebyshev polynomial for the entire
925 [0,1] y interval, we break the interval up into 100 equal 991 [0,1] y interval, we break the interval up into 100 equal
926 subintervals, with a switch/lookup table, and use much lower 992 subintervals, with a switch/lookup table, and use much lower
927 degree Chebyshev polynomials in each subinterval. This greatly 993 degree Chebyshev polynomials in each subinterval. This greatly
928 improves performance in my tests. 994 improves performance in my tests.
929 995
930 For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x), 996 For x < 0, we use the relationship erfcx(-x) = 2 exp(x^2) - erfc(x),
931 with the usual checks for overflow etcetera. 997 with the usual checks for overflow etcetera.
932 998
933 Performance-wise, it seems to be substantially faster than either 999 Performance-wise, it seems to be substantially faster than either
1355 { 1421 {
1356 if (x >= 0) { 1422 if (x >= 0) {
1357 if (x > 50) { // continued-fraction expansion is faster 1423 if (x > 50) { // continued-fraction expansion is faster
1358 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1424 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1359 if (x > 5e7) // 1-term expansion, important to avoid overflow 1425 if (x > 5e7) // 1-term expansion, important to avoid overflow
1360 return ispi / x; 1426 return ispi / x;
1361 /* 5-term expansion (rely on compiler for CSE), simplified from: 1427 /* 5-term expansion (rely on compiler for CSE), simplified from:
1362 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */ 1428 ispi / (x+0.5/(x+1/(x+1.5/(x+2/x)))) */
1363 return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75)); 1429 return ispi*((x*x) * (x*x+4.5) + 2) / (x * ((x*x) * (x*x+5) + 3.75));
1364 } 1430 }
1365 return erfcx_y100(400/(4+x)); 1431 return erfcx_y100(400/(4+x));
1366 } 1432 }
1367 else 1433 else
1368 return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x) 1434 return x < -26.7 ? HUGE_VAL : (x < -6.1 ? 2*exp(x*x)
1369 : 2*exp(x*x) - erfcx_y100(400/(4-x))); 1435 : 2*exp(x*x) - erfcx_y100(400/(4-x)));
1370 } 1436 }
1371 1437
1372 ///////////////////////////////////////////////////////////////////////// 1438 /////////////////////////////////////////////////////////////////////////
1373 /* Compute a scaled Dawson integral 1439 /* Compute a scaled Dawson integral
1374 FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi) 1440 FADDEEVA(w_im)(x) = 2*Dawson(x)/sqrt(pi)
1775 } 1841 }
1776 case 96: { 1842 case 96: {
1777 double t = 2*y100 - 193; 1843 double t = 2*y100 - 193;
1778 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; 1844 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;
1779 } 1845 }
1780 case 97: { 1846 case 97: case 98:
1781 double t = 2*y100 - 195; 1847 case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
1782 return 0.28920121009594899986e-1 + (-0.59271325915413781788e-2 + (0.28796136372768177423e-4 + (-0.30300382596279568642e-7 + (-0.97688275022802329749e-9 + (0.12179215701512592356e-10 - 0.93380988481777777779e-13 * t) * t) * t) * t) * t) * t; 1848 // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7 + 16/945 x^9)
1783 }
1784 case 98: {
1785 double t = 2*y100 - 197;
1786 return 0.17180782722617876655e-1 + (-0.58123419543161127769e-2 + (0.28591841095380959666e-4 + (-0.37642963496443667043e-7 + (-0.86055809047367300024e-9 + (0.11101709356762665578e-10 - 0.86272947493333333334e-13 * t) * t) * t) * t) * t) * t;
1787 }
1788 case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.010101...)
1789 // (2/sqrt(pi)) * (x - 2/3 x^3 + 4/15 x^5 - 8/105 x^7)
1790 double x2 = x*x; 1849 double x2 = x*x;
1791 return x * (1.1283791670955125739 1850 return x * (1.1283791670955125739
1792 - x2 * (0.75225277806367504925 1851 - x2 * (0.75225277806367504925
1793 - x2 * (0.30090111122547001970 1852 - x2 * (0.30090111122547001970
1794 - x2 * 0.085971746064420005629))); 1853 - x2 * (0.085971746064420005629
1854 - x2 * 0.016931216931216931217))));
1795 } 1855 }
1796 } 1856 }
1797 /* Since 0 <= y100 < 101, this is only reached if x is NaN, 1857 /* Since 0 <= y100 < 101, this is only reached if x is NaN,
1798 in which case we should return NaN. */ 1858 in which case we should return NaN. */
1799 return NaN; 1859 return NaN;
1803 { 1863 {
1804 if (x >= 0) { 1864 if (x >= 0) {
1805 if (x > 45) { // continued-fraction expansion is faster 1865 if (x > 45) { // continued-fraction expansion is faster
1806 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1866 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1807 if (x > 5e7) // 1-term expansion, important to avoid overflow 1867 if (x > 5e7) // 1-term expansion, important to avoid overflow
1808 return ispi / x; 1868 return ispi / x;
1809 /* 5-term expansion (rely on compiler for CSE), simplified from: 1869 /* 5-term expansion (rely on compiler for CSE), simplified from:
1810 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ 1870 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1811 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75)); 1871 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1812 } 1872 }
1813 return w_im_y100(100/(1+x), x); 1873 return w_im_y100(100/(1+x), x);
1814 } 1874 }
1815 else { // = -FADDEEVA(w_im)(-x) 1875 else { // = -FADDEEVA(w_im)(-x)
1816 if (x < -45) { // continued-fraction expansion is faster 1876 if (x < -45) { // continued-fraction expansion is faster
1817 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) 1877 const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)
1818 if (x < -5e7) // 1-term expansion, important to avoid overflow 1878 if (x < -5e7) // 1-term expansion, important to avoid overflow
1819 return ispi / x; 1879 return ispi / x;
1820 /* 5-term expansion (rely on compiler for CSE), simplified from: 1880 /* 5-term expansion (rely on compiler for CSE), simplified from:
1821 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */ 1881 ispi / (x-0.5/(x-1/(x-1.5/(x-2/x)))) */
1822 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75)); 1882 return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
1823 } 1883 }
1824 return -w_im_y100(100/(1-x), -x); 1884 return -w_im_y100(100/(1-x), -x);
1825 } 1885 }
1826 } 1886 }
1838 1898
1839 // compute relative error |b-a|/|a|, handling case of NaN and Inf, 1899 // compute relative error |b-a|/|a|, handling case of NaN and Inf,
1840 static double relerr(double a, double b) { 1900 static double relerr(double a, double b) {
1841 if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) { 1901 if (isnan(a) || isnan(b) || isinf(a) || isinf(b)) {
1842 if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) || 1902 if ((isnan(a) && !isnan(b)) || (!isnan(a) && isnan(b)) ||
1843 (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) || 1903 (isinf(a) && !isinf(b)) || (!isinf(a) && isinf(b)) ||
1844 (isinf(a) && isinf(b) && a*b < 0)) 1904 (isinf(a) && isinf(b) && a*b < 0))
1845 return Inf; // "infinite" error 1905 return Inf; // "infinite" error
1846 return 0; // matching infinity/nan results counted as zero error 1906 return 0; // matching infinity/nan results counted as zero error
1847 } 1907 }
1848 if (a == 0) 1908 if (a == 0)
1849 return b == 0 ? 0 : Inf; 1909 return b == 0 ? 0 : Inf;
1855 double errmax_all = 0; 1915 double errmax_all = 0;
1856 { 1916 {
1857 printf("############# w(z) tests #############\n"); 1917 printf("############# w(z) tests #############\n");
1858 #define NTST 57 // define instead of const for C compatibility 1918 #define NTST 57 // define instead of const for C compatibility
1859 cmplx z[NTST] = { 1919 cmplx z[NTST] = {
1860 CMPLX(624.2,-0.26123), 1920 C(624.2,-0.26123),
1861 CMPLX(-0.4,3.), 1921 C(-0.4,3.),
1862 CMPLX(0.6,2.), 1922 C(0.6,2.),
1863 CMPLX(-1.,1.), 1923 C(-1.,1.),
1864 CMPLX(-1.,-9.), 1924 C(-1.,-9.),
1865 CMPLX(-1.,9.), 1925 C(-1.,9.),
1866 CMPLX(-0.0000000234545,1.1234), 1926 C(-0.0000000234545,1.1234),
1867 CMPLX(-3.,5.1), 1927 C(-3.,5.1),
1868 CMPLX(-53,30.1), 1928 C(-53,30.1),
1869 CMPLX(0.0,0.12345), 1929 C(0.0,0.12345),
1870 CMPLX(11,1), 1930 C(11,1),
1871 CMPLX(-22,-2), 1931 C(-22,-2),
1872 CMPLX(9,-28), 1932 C(9,-28),
1873 CMPLX(21,-33), 1933 C(21,-33),
1874 CMPLX(1e5,1e5), 1934 C(1e5,1e5),
1875 CMPLX(1e14,1e14), 1935 C(1e14,1e14),
1876 CMPLX(-3001,-1000), 1936 C(-3001,-1000),
1877 CMPLX(1e160,-1e159), 1937 C(1e160,-1e159),
1878 CMPLX(-6.01,0.01), 1938 C(-6.01,0.01),
1879 CMPLX(-0.7,-0.7), 1939 C(-0.7,-0.7),
1880 CMPLX(2.611780000000000e+01, 4.540909610972489e+03), 1940 C(2.611780000000000e+01, 4.540909610972489e+03),
1881 CMPLX(0.8e7,0.3e7), 1941 C(0.8e7,0.3e7),
1882 CMPLX(-20,-19.8081), 1942 C(-20,-19.8081),
1883 CMPLX(1e-16,-1.1e-16), 1943 C(1e-16,-1.1e-16),
1884 CMPLX(2.3e-8,1.3e-8), 1944 C(2.3e-8,1.3e-8),
1885 CMPLX(6.3,-1e-13), 1945 C(6.3,-1e-13),
1886 CMPLX(6.3,1e-20), 1946 C(6.3,1e-20),
1887 CMPLX(1e-20,6.3), 1947 C(1e-20,6.3),
1888 CMPLX(1e-20,16.3), 1948 C(1e-20,16.3),
1889 CMPLX(9,1e-300), 1949 C(9,1e-300),
1890 CMPLX(6.01,0.11), 1950 C(6.01,0.11),
1891 CMPLX(8.01,1.01e-10), 1951 C(8.01,1.01e-10),
1892 CMPLX(28.01,1e-300), 1952 C(28.01,1e-300),
1893 CMPLX(10.01,1e-200), 1953 C(10.01,1e-200),
1894 CMPLX(10.01,-1e-200), 1954 C(10.01,-1e-200),
1895 CMPLX(10.01,0.99e-10), 1955 C(10.01,0.99e-10),
1896 CMPLX(10.01,-0.99e-10), 1956 C(10.01,-0.99e-10),
1897 CMPLX(1e-20,7.01), 1957 C(1e-20,7.01),
1898 CMPLX(-1,7.01), 1958 C(-1,7.01),
1899 CMPLX(5.99,7.01), 1959 C(5.99,7.01),
1900 CMPLX(1,0), 1960 C(1,0),
1901 CMPLX(55,0), 1961 C(55,0),
1902 CMPLX(-0.1,0), 1962 C(-0.1,0),
1903 CMPLX(1e-20,0), 1963 C(1e-20,0),
1904 CMPLX(0,5e-14), 1964 C(0,5e-14),
1905 CMPLX(0,51), 1965 C(0,51),
1906 CMPLX(Inf,0), 1966 C(Inf,0),
1907 CMPLX(-Inf,0), 1967 C(-Inf,0),
1908 CMPLX(0,Inf), 1968 C(0,Inf),
1909 CMPLX(0,-Inf), 1969 C(0,-Inf),
1910 CMPLX(Inf,Inf), 1970 C(Inf,Inf),
1911 CMPLX(Inf,-Inf), 1971 C(Inf,-Inf),
1912 CMPLX(NaN,NaN), 1972 C(NaN,NaN),
1913 CMPLX(NaN,0), 1973 C(NaN,0),
1914 CMPLX(0,NaN), 1974 C(0,NaN),
1915 CMPLX(NaN,Inf), 1975 C(NaN,Inf),
1916 CMPLX(Inf,NaN) 1976 C(Inf,NaN)
1917 }; 1977 };
1918 cmplx w[NTST] = { /* w(z), computed with WolframAlpha 1978 cmplx w[NTST] = { /* w(z), computed with WolframAlpha
1919 ... note that WolframAlpha is problematic 1979 ... note that WolframAlpha is problematic
1920 some of the above inputs, so I had to 1980 some of the above inputs, so I had to
1921 use the continued-fraction expansion 1981 use the continued-fraction expansion
1922 in WolframAlpha in some cases, or switch 1982 in WolframAlpha in some cases, or switch
1923 to Maple */ 1983 to Maple */
1924 CMPLX(-3.78270245518980507452677445620103199303131110e-7, 1984 C(-3.78270245518980507452677445620103199303131110e-7,
1925 0.000903861276433172057331093754199933411710053155), 1985 0.000903861276433172057331093754199933411710053155),
1926 CMPLX(0.1764906227004816847297495349730234591778719532788, 1986 C(0.1764906227004816847297495349730234591778719532788,
1927 -0.02146550539468457616788719893991501311573031095617), 1987 -0.02146550539468457616788719893991501311573031095617),
1928 CMPLX(0.2410250715772692146133539023007113781272362309451, 1988 C(0.2410250715772692146133539023007113781272362309451,
1929 0.06087579663428089745895459735240964093522265589350), 1989 0.06087579663428089745895459735240964093522265589350),
1930 CMPLX(0.30474420525691259245713884106959496013413834051768, 1990 C(0.30474420525691259245713884106959496013413834051768,
1931 -0.20821893820283162728743734725471561394145872072738), 1991 -0.20821893820283162728743734725471561394145872072738),
1932 CMPLX(7.317131068972378096865595229600561710140617977e34, 1992 C(7.317131068972378096865595229600561710140617977e34,
1933 8.321873499714402777186848353320412813066170427e34), 1993 8.321873499714402777186848353320412813066170427e34),
1934 CMPLX(0.0615698507236323685519612934241429530190806818395, 1994 C(0.0615698507236323685519612934241429530190806818395,
1935 -0.00676005783716575013073036218018565206070072304635), 1995 -0.00676005783716575013073036218018565206070072304635),
1936 CMPLX(0.3960793007699874918961319170187598400134746631, 1996 C(0.3960793007699874918961319170187598400134746631,
1937 -5.593152259116644920546186222529802777409274656e-9), 1997 -5.593152259116644920546186222529802777409274656e-9),
1938 CMPLX(0.08217199226739447943295069917990417630675021771804, 1998 C(0.08217199226739447943295069917990417630675021771804,
1939 -0.04701291087643609891018366143118110965272615832184), 1999 -0.04701291087643609891018366143118110965272615832184),
1940 CMPLX(0.00457246000350281640952328010227885008541748668738, 2000 C(0.00457246000350281640952328010227885008541748668738,
1941 -0.00804900791411691821818731763401840373998654987934), 2001 -0.00804900791411691821818731763401840373998654987934),
1942 CMPLX(0.8746342859608052666092782112565360755791467973338452, 2002 C(0.8746342859608052666092782112565360755791467973338452,
1943 0.), 2003 0.),
1944 CMPLX(0.00468190164965444174367477874864366058339647648741, 2004 C(0.00468190164965444174367477874864366058339647648741,
1945 0.0510735563901306197993676329845149741675029197050), 2005 0.0510735563901306197993676329845149741675029197050),
1946 CMPLX(-0.0023193175200187620902125853834909543869428763219, 2006 C(-0.0023193175200187620902125853834909543869428763219,
1947 -0.025460054739731556004902057663500272721780776336), 2007 -0.025460054739731556004902057663500272721780776336),
1948 CMPLX(9.11463368405637174660562096516414499772662584e304, 2008 C(9.11463368405637174660562096516414499772662584e304,
1949 3.97101807145263333769664875189354358563218932e305), 2009 3.97101807145263333769664875189354358563218932e305),
1950 CMPLX(-4.4927207857715598976165541011143706155432296e281, 2010 C(-4.4927207857715598976165541011143706155432296e281,
1951 -2.8019591213423077494444700357168707775769028e281), 2011 -2.8019591213423077494444700357168707775769028e281),
1952 CMPLX(2.820947917809305132678577516325951485807107151e-6, 2012 C(2.820947917809305132678577516325951485807107151e-6,
1953 2.820947917668257736791638444590253942253354058e-6), 2013 2.820947917668257736791638444590253942253354058e-6),
1954 CMPLX(2.82094791773878143474039725787438662716372268e-15, 2014 C(2.82094791773878143474039725787438662716372268e-15,
1955 2.82094791773878143474039725773333923127678361e-15), 2015 2.82094791773878143474039725773333923127678361e-15),
1956 CMPLX(-0.0000563851289696244350147899376081488003110150498, 2016 C(-0.0000563851289696244350147899376081488003110150498,
1957 -0.000169211755126812174631861529808288295454992688), 2017 -0.000169211755126812174631861529808288295454992688),
1958 CMPLX(-5.586035480670854326218608431294778077663867e-162, 2018 C(-5.586035480670854326218608431294778077663867e-162,
1959 5.586035480670854326218608431294778077663867e-161), 2019 5.586035480670854326218608431294778077663867e-161),
1960 CMPLX(0.00016318325137140451888255634399123461580248456, 2020 C(0.00016318325137140451888255634399123461580248456,
1961 -0.095232456573009287370728788146686162555021209999), 2021 -0.095232456573009287370728788146686162555021209999),
1962 CMPLX(0.69504753678406939989115375989939096800793577783885, 2022 C(0.69504753678406939989115375989939096800793577783885,
1963 -1.8916411171103639136680830887017670616339912024317), 2023 -1.8916411171103639136680830887017670616339912024317),
1964 CMPLX(0.0001242418269653279656612334210746733213167234822, 2024 C(0.0001242418269653279656612334210746733213167234822,
1965 7.145975826320186888508563111992099992116786763e-7), 2025 7.145975826320186888508563111992099992116786763e-7),
1966 CMPLX(2.318587329648353318615800865959225429377529825e-8, 2026 C(2.318587329648353318615800865959225429377529825e-8,
1967 6.182899545728857485721417893323317843200933380e-8), 2027 6.182899545728857485721417893323317843200933380e-8),
1968 CMPLX(-0.0133426877243506022053521927604277115767311800303, 2028 C(-0.0133426877243506022053521927604277115767311800303,
1969 -0.0148087097143220769493341484176979826888871576145), 2029 -0.0148087097143220769493341484176979826888871576145),
1970 CMPLX(1.00000000000000012412170838050638522857747934, 2030 C(1.00000000000000012412170838050638522857747934,
1971 1.12837916709551279389615890312156495593616433e-16), 2031 1.12837916709551279389615890312156495593616433e-16),
1972 CMPLX(0.9999999853310704677583504063775310832036830015, 2032 C(0.9999999853310704677583504063775310832036830015,
1973 2.595272024519678881897196435157270184030360773e-8), 2033 2.595272024519678881897196435157270184030360773e-8),
1974 CMPLX(-1.4731421795638279504242963027196663601154624e-15, 2034 C(-1.4731421795638279504242963027196663601154624e-15,
1975 0.090727659684127365236479098488823462473074709), 2035 0.090727659684127365236479098488823462473074709),
1976 CMPLX(5.79246077884410284575834156425396800754409308e-18, 2036 C(5.79246077884410284575834156425396800754409308e-18,
1977 0.0907276596841273652364790985059772809093822374), 2037 0.0907276596841273652364790985059772809093822374),
1978 CMPLX(0.0884658993528521953466533278764830881245144368, 2038 C(0.0884658993528521953466533278764830881245144368,
1979 1.37088352495749125283269718778582613192166760e-22), 2039 1.37088352495749125283269718778582613192166760e-22),
1980 CMPLX(0.0345480845419190424370085249304184266813447878, 2040 C(0.0345480845419190424370085249304184266813447878,
1981 2.11161102895179044968099038990446187626075258e-23), 2041 2.11161102895179044968099038990446187626075258e-23),
1982 CMPLX(6.63967719958073440070225527042829242391918213e-36, 2042 C(6.63967719958073440070225527042829242391918213e-36,
1983 0.0630820900592582863713653132559743161572639353), 2043 0.0630820900592582863713653132559743161572639353),
1984 CMPLX(0.00179435233208702644891092397579091030658500743634, 2044 C(0.00179435233208702644891092397579091030658500743634,
1985 0.0951983814805270647939647438459699953990788064762), 2045 0.0951983814805270647939647438459699953990788064762),
1986 CMPLX(9.09760377102097999924241322094863528771095448e-13, 2046 C(9.09760377102097999924241322094863528771095448e-13,
1987 0.0709979210725138550986782242355007611074966717), 2047 0.0709979210725138550986782242355007611074966717),
1988 CMPLX(7.2049510279742166460047102593255688682910274423e-304, 2048 C(7.2049510279742166460047102593255688682910274423e-304,
1989 0.0201552956479526953866611812593266285000876784321), 2049 0.0201552956479526953866611812593266285000876784321),
1990 CMPLX(3.04543604652250734193622967873276113872279682e-44, 2050 C(3.04543604652250734193622967873276113872279682e-44,
1991 0.0566481651760675042930042117726713294607499165), 2051 0.0566481651760675042930042117726713294607499165),
1992 CMPLX(3.04543604652250734193622967873276113872279682e-44, 2052 C(3.04543604652250734193622967873276113872279682e-44,
1993 0.0566481651760675042930042117726713294607499165), 2053 0.0566481651760675042930042117726713294607499165),
1994 CMPLX(0.5659928732065273429286988428080855057102069081e-12, 2054 C(0.5659928732065273429286988428080855057102069081e-12,
1995 0.056648165176067504292998527162143030538756683302), 2055 0.056648165176067504292998527162143030538756683302),
1996 CMPLX(-0.56599287320652734292869884280802459698927645e-12, 2056 C(-0.56599287320652734292869884280802459698927645e-12,
1997 0.0566481651760675042929985271621430305387566833029), 2057 0.0566481651760675042929985271621430305387566833029),
1998 CMPLX(0.0796884251721652215687859778119964009569455462, 2058 C(0.0796884251721652215687859778119964009569455462,
1999 1.11474461817561675017794941973556302717225126e-22), 2059 1.11474461817561675017794941973556302717225126e-22),
2000 CMPLX(0.07817195821247357458545539935996687005781943386550, 2060 C(0.07817195821247357458545539935996687005781943386550,
2001 -0.01093913670103576690766705513142246633056714279654), 2061 -0.01093913670103576690766705513142246633056714279654),
2002 CMPLX(0.04670032980990449912809326141164730850466208439937, 2062 C(0.04670032980990449912809326141164730850466208439937,
2003 0.03944038961933534137558064191650437353429669886545), 2063 0.03944038961933534137558064191650437353429669886545),
2004 CMPLX(0.36787944117144232159552377016146086744581113103176, 2064 C(0.36787944117144232159552377016146086744581113103176,
2005 0.60715770584139372911503823580074492116122092866515), 2065 0.60715770584139372911503823580074492116122092866515),
2006 CMPLX(0, 2066 C(0,
2007 0.010259688805536830986089913987516716056946786526145), 2067 0.010259688805536830986089913987516716056946786526145),
2008 CMPLX(0.99004983374916805357390597718003655777207908125383, 2068 C(0.99004983374916805357390597718003655777207908125383,
2009 -0.11208866436449538036721343053869621153527769495574), 2069 -0.11208866436449538036721343053869621153527769495574),
2010 CMPLX(0.99999999999999999999999999999999999999990000, 2070 C(0.99999999999999999999999999999999999999990000,
2011 1.12837916709551257389615890312154517168802603e-20), 2071 1.12837916709551257389615890312154517168802603e-20),
2012 CMPLX(0.999999999999943581041645226871305192054749891144158, 2072 C(0.999999999999943581041645226871305192054749891144158,
2013 0), 2073 0),
2014 CMPLX(0.0110604154853277201542582159216317923453996211744250, 2074 C(0.0110604154853277201542582159216317923453996211744250,
2015 0), 2075 0),
2016 CMPLX(0,0), 2076 C(0,0),
2017 CMPLX(0,0), 2077 C(0,0),
2018 CMPLX(0,0), 2078 C(0,0),
2019 CMPLX(Inf,0), 2079 C(Inf,0),
2020 CMPLX(0,0), 2080 C(0,0),
2021 CMPLX(NaN,NaN), 2081 C(NaN,NaN),
2022 CMPLX(NaN,NaN), 2082 C(NaN,NaN),
2023 CMPLX(NaN,NaN), 2083 C(NaN,NaN),
2024 CMPLX(NaN,0), 2084 C(NaN,0),
2025 CMPLX(NaN,NaN), 2085 C(NaN,NaN),
2026 CMPLX(NaN,NaN) 2086 C(NaN,NaN)
2027 }; 2087 };
2028 double errmax = 0; 2088 double errmax = 0;
2029 for (int i = 0; i < NTST; ++i) { 2089 for (int i = 0; i < NTST; ++i) {
2030 cmplx fw = FADDEEVA(w)(z[i],0.); 2090 cmplx fw = FADDEEVA(w)(z[i],0.);
2031 double re_err = relerr(creal(w[i]), creal(fw)); 2091 double re_err = relerr(creal(w[i]), creal(fw));
2032 double im_err = relerr(cimag(w[i]), cimag(fw)); 2092 double im_err = relerr(cimag(w[i]), cimag(fw));
2033 printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", 2093 printf("w(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n",
2034 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), 2094 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]),
2035 re_err, im_err); 2095 re_err, im_err);
2036 if (re_err > errmax) errmax = re_err; 2096 if (re_err > errmax) errmax = re_err;
2037 if (im_err > errmax) errmax = im_err; 2097 if (im_err > errmax) errmax = im_err;
2038 } 2098 }
2039 if (errmax > 1e-13) { 2099 if (errmax > 1e-13) {
2040 printf("FAILURE -- relative error %g too large!\n", errmax); 2100 printf("FAILURE -- relative error %g too large!\n", errmax);
2043 printf("SUCCESS (max relative error = %g)\n", errmax); 2103 printf("SUCCESS (max relative error = %g)\n", errmax);
2044 if (errmax > errmax_all) errmax_all = errmax; 2104 if (errmax > errmax_all) errmax_all = errmax;
2045 } 2105 }
2046 { 2106 {
2047 #undef NTST 2107 #undef NTST
2048 #define NTST 33 // define instead of const for C compatibility 2108 #define NTST 41 // define instead of const for C compatibility
2049 cmplx z[NTST] = { 2109 cmplx z[NTST] = {
2050 CMPLX(1,2), 2110 C(1,2),
2051 CMPLX(-1,2), 2111 C(-1,2),
2052 CMPLX(1,-2), 2112 C(1,-2),
2053 CMPLX(-1,-2), 2113 C(-1,-2),
2054 CMPLX(9,-28), 2114 C(9,-28),
2055 CMPLX(21,-33), 2115 C(21,-33),
2056 CMPLX(1e3,1e3), 2116 C(1e3,1e3),
2057 CMPLX(-3001,-1000), 2117 C(-3001,-1000),
2058 CMPLX(1e160,-1e159), 2118 C(1e160,-1e159),
2059 CMPLX(5.1e-3, 1e-8), 2119 C(5.1e-3, 1e-8),
2060 CMPLX(-4.9e-3, 4.95e-3), 2120 C(-4.9e-3, 4.95e-3),
2061 CMPLX(4.9e-3, 0.5), 2121 C(4.9e-3, 0.5),
2062 CMPLX(4.9e-4, -0.5e1), 2122 C(4.9e-4, -0.5e1),
2063 CMPLX(-4.9e-5, -0.5e2), 2123 C(-4.9e-5, -0.5e2),
2064 CMPLX(5.1e-3, 0.5), 2124 C(5.1e-3, 0.5),
2065 CMPLX(5.1e-4, -0.5e1), 2125 C(5.1e-4, -0.5e1),
2066 CMPLX(-5.1e-5, -0.5e2), 2126 C(-5.1e-5, -0.5e2),
2067 CMPLX(1e-6,2e-6), 2127 C(1e-6,2e-6),
2068 CMPLX(0,2e-6), 2128 C(0,2e-6),
2069 CMPLX(0,2), 2129 C(0,2),
2070 CMPLX(0,20), 2130 C(0,20),
2071 CMPLX(0,200), 2131 C(0,200),
2072 CMPLX(Inf,0), 2132 C(Inf,0),
2073 CMPLX(-Inf,0), 2133 C(-Inf,0),
2074 CMPLX(0,Inf), 2134 C(0,Inf),
2075 CMPLX(0,-Inf), 2135 C(0,-Inf),
2076 CMPLX(Inf,Inf), 2136 C(Inf,Inf),
2077 CMPLX(Inf,-Inf), 2137 C(Inf,-Inf),
2078 CMPLX(NaN,NaN), 2138 C(NaN,NaN),
2079 CMPLX(NaN,0), 2139 C(NaN,0),
2080 CMPLX(0,NaN), 2140 C(0,NaN),
2081 CMPLX(NaN,Inf), 2141 C(NaN,Inf),
2082 CMPLX(Inf,NaN) 2142 C(Inf,NaN),
2143 C(1e-3,NaN),
2144 C(7e-2,7e-2),
2145 C(7e-2,-7e-4),
2146 C(-9e-2,7e-4),
2147 C(-9e-2,9e-2),
2148 C(-7e-4,9e-2),
2149 C(7e-2,0.9e-2),
2150 C(7e-2,1.1e-2)
2083 }; 2151 };
2084 cmplx w[NTST] = { // erf(z[i]), evaluated with Maple 2152 cmplx w[NTST] = { // erf(z[i]), evaluated with Maple
2085 CMPLX(-0.5366435657785650339917955593141927494421, 2153 C(-0.5366435657785650339917955593141927494421,
2086 -5.049143703447034669543036958614140565553), 2154 -5.049143703447034669543036958614140565553),
2087 CMPLX(0.5366435657785650339917955593141927494421, 2155 C(0.5366435657785650339917955593141927494421,
2088 -5.049143703447034669543036958614140565553), 2156 -5.049143703447034669543036958614140565553),
2089 CMPLX(-0.5366435657785650339917955593141927494421, 2157 C(-0.5366435657785650339917955593141927494421,
2090 5.049143703447034669543036958614140565553), 2158 5.049143703447034669543036958614140565553),
2091 CMPLX(0.5366435657785650339917955593141927494421, 2159 C(0.5366435657785650339917955593141927494421,
2092 5.049143703447034669543036958614140565553), 2160 5.049143703447034669543036958614140565553),
2093 CMPLX(0.3359473673830576996788000505817956637777e304, 2161 C(0.3359473673830576996788000505817956637777e304,
2094 -0.1999896139679880888755589794455069208455e304), 2162 -0.1999896139679880888755589794455069208455e304),
2095 CMPLX(0.3584459971462946066523939204836760283645e278, 2163 C(0.3584459971462946066523939204836760283645e278,
2096 0.3818954885257184373734213077678011282505e280), 2164 0.3818954885257184373734213077678011282505e280),
2097 CMPLX(0.9996020422657148639102150147542224526887, 2165 C(0.9996020422657148639102150147542224526887,
2098 0.00002801044116908227889681753993542916894856), 2166 0.00002801044116908227889681753993542916894856),
2099 CMPLX(-1, 0), 2167 C(-1, 0),
2100 CMPLX(1, 0), 2168 C(1, 0),
2101 CMPLX(0.005754683859034800134412990541076554934877, 2169 C(0.005754683859034800134412990541076554934877,
2102 0.1128349818335058741511924929801267822634e-7), 2170 0.1128349818335058741511924929801267822634e-7),
2103 CMPLX(-0.005529149142341821193633460286828381876955, 2171 C(-0.005529149142341821193633460286828381876955,
2104 0.005585388387864706679609092447916333443570), 2172 0.005585388387864706679609092447916333443570),
2105 CMPLX(0.007099365669981359632319829148438283865814, 2173 C(0.007099365669981359632319829148438283865814,
2106 0.6149347012854211635026981277569074001219), 2174 0.6149347012854211635026981277569074001219),
2107 CMPLX(0.3981176338702323417718189922039863062440e8, 2175 C(0.3981176338702323417718189922039863062440e8,
2108 -0.8298176341665249121085423917575122140650e10), 2176 -0.8298176341665249121085423917575122140650e10),
2109 CMPLX(-Inf, 2177 C(-Inf,
2110 -Inf), 2178 -Inf),
2111 CMPLX(0.007389128308257135427153919483147229573895, 2179 C(0.007389128308257135427153919483147229573895,
2112 0.6149332524601658796226417164791221815139), 2180 0.6149332524601658796226417164791221815139),
2113 CMPLX(0.4143671923267934479245651547534414976991e8, 2181 C(0.4143671923267934479245651547534414976991e8,
2114 -0.8298168216818314211557046346850921446950e10), 2182 -0.8298168216818314211557046346850921446950e10),
2115 CMPLX(-Inf, 2183 C(-Inf,
2116 -Inf), 2184 -Inf),
2117 CMPLX(0.1128379167099649964175513742247082845155e-5, 2185 C(0.1128379167099649964175513742247082845155e-5,
2118 0.2256758334191777400570377193451519478895e-5), 2186 0.2256758334191777400570377193451519478895e-5),
2119 CMPLX(0, 2187 C(0,
2120 0.2256758334194034158904576117253481476197e-5), 2188 0.2256758334194034158904576117253481476197e-5),
2121 CMPLX(0, 2189 C(0,
2122 18.56480241457555259870429191324101719886), 2190 18.56480241457555259870429191324101719886),
2123 CMPLX(0, 2191 C(0,
2124 0.1474797539628786202447733153131835124599e173), 2192 0.1474797539628786202447733153131835124599e173),
2125 CMPLX(0, 2193 C(0,
2126 Inf), 2194 Inf),
2127 CMPLX(1,0), 2195 C(1,0),
2128 CMPLX(-1,0), 2196 C(-1,0),
2129 CMPLX(0,Inf), 2197 C(0,Inf),
2130 CMPLX(0,-Inf), 2198 C(0,-Inf),
2131 CMPLX(NaN,NaN), 2199 C(NaN,NaN),
2132 CMPLX(NaN,NaN), 2200 C(NaN,NaN),
2133 CMPLX(NaN,NaN), 2201 C(NaN,NaN),
2134 CMPLX(NaN,0), 2202 C(NaN,0),
2135 CMPLX(0,NaN), 2203 C(0,NaN),
2136 CMPLX(NaN,NaN), 2204 C(NaN,NaN),
2137 CMPLX(NaN,NaN) 2205 C(NaN,NaN),
2206 C(NaN,NaN),
2207 C(0.07924380404615782687930591956705225541145,
2208 0.07872776218046681145537914954027729115247),
2209 C(0.07885775828512276968931773651224684454495,
2210 -0.0007860046704118224342390725280161272277506),
2211 C(-0.1012806432747198859687963080684978759881,
2212 0.0007834934747022035607566216654982820299469),
2213 C(-0.1020998418798097910247132140051062512527,
2214 0.1010030778892310851309082083238896270340),
2215 C(-0.0007962891763147907785684591823889484764272,
2216 0.1018289385936278171741809237435404896152),
2217 C(0.07886408666470478681566329888615410479530,
2218 0.01010604288780868961492224347707949372245),
2219 C(0.07886723099940260286824654364807981336591,
2220 0.01235199327873258197931147306290916629654)
2138 }; 2221 };
2139 #define TST(f,isc) \ 2222 #define TST(f,isc) \
2140 printf("############# " #f "(z) tests #############\n"); \ 2223 printf("############# " #f "(z) tests #############\n"); \
2141 double errmax = 0; \ 2224 double errmax = 0; \
2142 for (int i = 0; i < NTST; ++i) { \ 2225 for (int i = 0; i < NTST; ++i) { \
2143 cmplx fw = FADDEEVA(f)(z[i],0.); \ 2226 cmplx fw = FADDEEVA(f)(z[i],0.); \
2144 double re_err = relerr(creal(w[i]), creal(fw)); \ 2227 double re_err = relerr(creal(w[i]), creal(fw)); \
2145 double im_err = relerr(cimag(w[i]), cimag(fw)); \ 2228 double im_err = relerr(cimag(w[i]), cimag(fw)); \
2146 printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \ 2229 printf(#f "(%g%+gi) = %g%+gi (vs. %g%+gi), re/im rel. err. = %0.2g/%0.2g)\n", \
2147 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \ 2230 creal(z[i]),cimag(z[i]), creal(fw),cimag(fw), creal(w[i]),cimag(w[i]), \
2148 re_err, im_err); \ 2231 re_err, im_err); \
2149 if (re_err > errmax) errmax = re_err; \ 2232 if (re_err > errmax) errmax = re_err; \
2150 if (im_err > errmax) errmax = im_err; \ 2233 if (im_err > errmax) errmax = im_err; \
2151 } \ 2234 } \
2152 if (errmax > 1e-13) { \ 2235 if (errmax > 1e-13) { \
2153 printf("FAILURE -- relative error %g too large!\n", errmax); \ 2236 printf("FAILURE -- relative error %g too large!\n", errmax); \
2154 return 1; \ 2237 return 1; \
2155 } \ 2238 } \
2156 printf("Checking " #f "(x) special case...\n"); \ 2239 printf("Checking " #f "(x) special case...\n"); \
2157 for (int i = 0; i < 10000; ++i) { \ 2240 for (int i = 0; i < 10000; ++i) { \
2158 double x = pow(10., -300. + i * 600. / (10000 - 1)); \ 2241 double x = pow(10., -300. + i * 600. / (10000 - 1)); \
2159 double re_err = relerr(FADDEEVA_RE(f)(x), \ 2242 double re_err = relerr(FADDEEVA_RE(f)(x), \
2160 creal(FADDEEVA(f)(CMPLX(x,x*isc),0.))); \ 2243 creal(FADDEEVA(f)(C(x,x*isc),0.))); \
2161 if (re_err > errmax) errmax = re_err; \ 2244 if (re_err > errmax) errmax = re_err; \
2162 re_err = relerr(FADDEEVA_RE(f)(-x), \ 2245 re_err = relerr(FADDEEVA_RE(f)(-x), \
2163 creal(FADDEEVA(f)(CMPLX(-x,x*isc),0.))); \ 2246 creal(FADDEEVA(f)(C(-x,x*isc),0.))); \
2164 if (re_err > errmax) errmax = re_err; \ 2247 if (re_err > errmax) errmax = re_err; \
2165 } \ 2248 } \
2166 { \ 2249 { \
2167 double re_err = relerr(FADDEEVA_RE(f)(Inf), \ 2250 double re_err = relerr(FADDEEVA_RE(f)(Inf), \
2168 creal(FADDEEVA(f)(CMPLX(Inf,0.),0.))); \ 2251 creal(FADDEEVA(f)(C(Inf,0.),0.))); \
2169 if (re_err > errmax) errmax = re_err; \ 2252 if (re_err > errmax) errmax = re_err; \
2170 re_err = relerr(FADDEEVA_RE(f)(-Inf), \ 2253 re_err = relerr(FADDEEVA_RE(f)(-Inf), \
2171 creal(FADDEEVA(f)(CMPLX(-Inf,0.),0.))); \ 2254 creal(FADDEEVA(f)(C(-Inf,0.),0.))); \
2172 if (re_err > errmax) errmax = re_err; \ 2255 if (re_err > errmax) errmax = re_err; \
2173 re_err = relerr(FADDEEVA_RE(f)(NaN), \ 2256 re_err = relerr(FADDEEVA_RE(f)(NaN), \
2174 creal(FADDEEVA(f)(CMPLX(NaN,0.),0.))); \ 2257 creal(FADDEEVA(f)(C(NaN,0.),0.))); \
2175 if (re_err > errmax) errmax = re_err; \ 2258 if (re_err > errmax) errmax = re_err; \
2176 } \ 2259 } \
2177 if (errmax > 1e-13) { \ 2260 if (errmax > 1e-13) { \
2178 printf("FAILURE -- relative error %g too large!\n", errmax); \ 2261 printf("FAILURE -- relative error %g too large!\n", errmax); \
2179 return 1; \ 2262 return 1; \
2180 } \ 2263 } \
2181 printf("SUCCESS (max relative error = %g)\n", errmax); \ 2264 printf("SUCCESS (max relative error = %g)\n", errmax); \
2182 if (errmax > errmax_all) errmax_all = errmax 2265 if (errmax > errmax_all) errmax_all = errmax
2183 2266
2184 TST(erf, 1e-20); 2267 TST(erf, 1e-20);
2185 } 2268 }
2186 { 2269 {
2187 // since erfi just calls through to erf, just one test should 2270 // since erfi just calls through to erf, just one test should
2188 // be sufficient to make sure I didn't screw up the signs or something 2271 // be sufficient to make sure I didn't screw up the signs or something
2189 #undef NTST 2272 #undef NTST
2190 #define NTST 1 // define instead of const for C compatibility 2273 #define NTST 1 // define instead of const for C compatibility
2191 cmplx z[NTST] = { CMPLX(1.234,0.5678) }; 2274 cmplx z[NTST] = { C(1.234,0.5678) };
2192 cmplx w[NTST] = { // erfi(z[i]), computed with Maple 2275 cmplx w[NTST] = { // erfi(z[i]), computed with Maple
2193 CMPLX(1.081032284405373149432716643834106923212, 2276 C(1.081032284405373149432716643834106923212,
2194 1.926775520840916645838949402886591180834) 2277 1.926775520840916645838949402886591180834)
2195 }; 2278 };
2196 TST(erfi, 0); 2279 TST(erfi, 0);
2197 } 2280 }
2198 { 2281 {
2199 // since erfcx just calls through to w, just one test should 2282 // since erfcx just calls through to w, just one test should
2200 // be sufficient to make sure I didn't screw up the signs or something 2283 // be sufficient to make sure I didn't screw up the signs or something
2201 #undef NTST 2284 #undef NTST
2202 #define NTST 1 // define instead of const for C compatibility 2285 #define NTST 1 // define instead of const for C compatibility
2203 cmplx z[NTST] = { CMPLX(1.234,0.5678) }; 2286 cmplx z[NTST] = { C(1.234,0.5678) };
2204 cmplx w[NTST] = { // erfcx(z[i]), computed with Maple 2287 cmplx w[NTST] = { // erfcx(z[i]), computed with Maple
2205 CMPLX(0.3382187479799972294747793561190487832579, 2288 C(0.3382187479799972294747793561190487832579,
2206 -0.1116077470811648467464927471872945833154) 2289 -0.1116077470811648467464927471872945833154)
2207 }; 2290 };
2208 TST(erfcx, 0); 2291 TST(erfcx, 0);
2209 } 2292 }
2210 { 2293 {
2211 #undef NTST 2294 #undef NTST
2212 #define NTST 30 // define instead of const for C compatibility 2295 #define NTST 30 // define instead of const for C compatibility
2213 cmplx z[NTST] = { 2296 cmplx z[NTST] = {
2214 CMPLX(1,2), 2297 C(1,2),
2215 CMPLX(-1,2), 2298 C(-1,2),
2216 CMPLX(1,-2), 2299 C(1,-2),
2217 CMPLX(-1,-2), 2300 C(-1,-2),
2218 CMPLX(9,-28), 2301 C(9,-28),
2219 CMPLX(21,-33), 2302 C(21,-33),
2220 CMPLX(1e3,1e3), 2303 C(1e3,1e3),
2221 CMPLX(-3001,-1000), 2304 C(-3001,-1000),
2222 CMPLX(1e160,-1e159), 2305 C(1e160,-1e159),
2223 CMPLX(5.1e-3, 1e-8), 2306 C(5.1e-3, 1e-8),
2224 CMPLX(0,2e-6), 2307 C(0,2e-6),
2225 CMPLX(0,2), 2308 C(0,2),
2226 CMPLX(0,20), 2309 C(0,20),
2227 CMPLX(0,200), 2310 C(0,200),
2228 CMPLX(2e-6,0), 2311 C(2e-6,0),
2229 CMPLX(2,0), 2312 C(2,0),
2230 CMPLX(20,0), 2313 C(20,0),
2231 CMPLX(200,0), 2314 C(200,0),
2232 CMPLX(Inf,0), 2315 C(Inf,0),
2233 CMPLX(-Inf,0), 2316 C(-Inf,0),
2234 CMPLX(0,Inf), 2317 C(0,Inf),
2235 CMPLX(0,-Inf), 2318 C(0,-Inf),
2236 CMPLX(Inf,Inf), 2319 C(Inf,Inf),
2237 CMPLX(Inf,-Inf), 2320 C(Inf,-Inf),
2238 CMPLX(NaN,NaN), 2321 C(NaN,NaN),
2239 CMPLX(NaN,0), 2322 C(NaN,0),
2240 CMPLX(0,NaN), 2323 C(0,NaN),
2241 CMPLX(NaN,Inf), 2324 C(NaN,Inf),
2242 CMPLX(Inf,NaN), 2325 C(Inf,NaN),
2243 CMPLX(88,0) 2326 C(88,0)
2244 }; 2327 };
2245 cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple 2328 cmplx w[NTST] = { // erfc(z[i]), evaluated with Maple
2246 CMPLX(1.536643565778565033991795559314192749442, 2329 C(1.536643565778565033991795559314192749442,
2247 5.049143703447034669543036958614140565553), 2330 5.049143703447034669543036958614140565553),
2248 CMPLX(0.4633564342214349660082044406858072505579, 2331 C(0.4633564342214349660082044406858072505579,
2249 5.049143703447034669543036958614140565553), 2332 5.049143703447034669543036958614140565553),
2250 CMPLX(1.536643565778565033991795559314192749442, 2333 C(1.536643565778565033991795559314192749442,
2251 -5.049143703447034669543036958614140565553), 2334 -5.049143703447034669543036958614140565553),
2252 CMPLX(0.4633564342214349660082044406858072505579, 2335 C(0.4633564342214349660082044406858072505579,
2253 -5.049143703447034669543036958614140565553), 2336 -5.049143703447034669543036958614140565553),
2254 CMPLX(-0.3359473673830576996788000505817956637777e304, 2337 C(-0.3359473673830576996788000505817956637777e304,
2255 0.1999896139679880888755589794455069208455e304), 2338 0.1999896139679880888755589794455069208455e304),
2256 CMPLX(-0.3584459971462946066523939204836760283645e278, 2339 C(-0.3584459971462946066523939204836760283645e278,
2257 -0.3818954885257184373734213077678011282505e280), 2340 -0.3818954885257184373734213077678011282505e280),
2258 CMPLX(0.0003979577342851360897849852457775473112748, 2341 C(0.0003979577342851360897849852457775473112748,
2259 -0.00002801044116908227889681753993542916894856), 2342 -0.00002801044116908227889681753993542916894856),
2260 CMPLX(2, 0), 2343 C(2, 0),
2261 CMPLX(0, 0), 2344 C(0, 0),
2262 CMPLX(0.9942453161409651998655870094589234450651, 2345 C(0.9942453161409651998655870094589234450651,
2263 -0.1128349818335058741511924929801267822634e-7), 2346 -0.1128349818335058741511924929801267822634e-7),
2264 CMPLX(1, 2347 C(1,
2265 -0.2256758334194034158904576117253481476197e-5), 2348 -0.2256758334194034158904576117253481476197e-5),
2266 CMPLX(1, 2349 C(1,
2267 -18.56480241457555259870429191324101719886), 2350 -18.56480241457555259870429191324101719886),
2268 CMPLX(1, 2351 C(1,
2269 -0.1474797539628786202447733153131835124599e173), 2352 -0.1474797539628786202447733153131835124599e173),
2270 CMPLX(1, -Inf), 2353 C(1, -Inf),
2271 CMPLX(0.9999977432416658119838633199332831406314, 2354 C(0.9999977432416658119838633199332831406314,
2272 0), 2355 0),
2273 CMPLX(0.004677734981047265837930743632747071389108, 2356 C(0.004677734981047265837930743632747071389108,
2274 0), 2357 0),
2275 CMPLX(0.5395865611607900928934999167905345604088e-175, 2358 C(0.5395865611607900928934999167905345604088e-175,
2276 0), 2359 0),
2277 CMPLX(0, 0), 2360 C(0, 0),
2278 CMPLX(0, 0), 2361 C(0, 0),
2279 CMPLX(2, 0), 2362 C(2, 0),
2280 CMPLX(1, -Inf), 2363 C(1, -Inf),
2281 CMPLX(1, Inf), 2364 C(1, Inf),
2282 CMPLX(NaN, NaN), 2365 C(NaN, NaN),
2283 CMPLX(NaN, NaN), 2366 C(NaN, NaN),
2284 CMPLX(NaN, NaN), 2367 C(NaN, NaN),
2285 CMPLX(NaN, 0), 2368 C(NaN, 0),
2286 CMPLX(1, NaN), 2369 C(1, NaN),
2287 CMPLX(NaN, NaN), 2370 C(NaN, NaN),
2288 CMPLX(NaN, NaN), 2371 C(NaN, NaN),
2289 CMPLX(0,0) 2372 C(0,0)
2290 }; 2373 };
2291 TST(erfc, 1e-20); 2374 TST(erfc, 1e-20);
2292 } 2375 }
2293 { 2376 {
2294 #undef NTST 2377 #undef NTST
2295 #define NTST 47 // define instead of const for C compatibility 2378 #define NTST 48 // define instead of const for C compatibility
2296 cmplx z[NTST] = { 2379 cmplx z[NTST] = {
2297 CMPLX(2,1), 2380 C(2,1),
2298 CMPLX(-2,1), 2381 C(-2,1),
2299 CMPLX(2,-1), 2382 C(2,-1),
2300 CMPLX(-2,-1), 2383 C(-2,-1),
2301 CMPLX(-28,9), 2384 C(-28,9),
2302 CMPLX(33,-21), 2385 C(33,-21),
2303 CMPLX(1e3,1e3), 2386 C(1e3,1e3),
2304 CMPLX(-1000,-3001), 2387 C(-1000,-3001),
2305 CMPLX(1e-8, 5.1e-3), 2388 C(1e-8, 5.1e-3),
2306 CMPLX(4.95e-3, -4.9e-3), 2389 C(4.95e-3, -4.9e-3),
2307 CMPLX(0.5, 4.9e-3), 2390 C(5.1e-3, 5.1e-3),
2308 CMPLX(-0.5e1, 4.9e-4), 2391 C(0.5, 4.9e-3),
2309 CMPLX(-0.5e2, -4.9e-5), 2392 C(-0.5e1, 4.9e-4),
2310 CMPLX(0.5e3, 4.9e-6), 2393 C(-0.5e2, -4.9e-5),
2311 CMPLX(0.5, 5.1e-3), 2394 C(0.5e3, 4.9e-6),
2312 CMPLX(-0.5e1, 5.1e-4), 2395 C(0.5, 5.1e-3),
2313 CMPLX(-0.5e2, -5.1e-5), 2396 C(-0.5e1, 5.1e-4),
2314 CMPLX(1e-6,2e-6), 2397 C(-0.5e2, -5.1e-5),
2315 CMPLX(2e-6,0), 2398 C(1e-6,2e-6),
2316 CMPLX(2,0), 2399 C(2e-6,0),
2317 CMPLX(20,0), 2400 C(2,0),
2318 CMPLX(200,0), 2401 C(20,0),
2319 CMPLX(0,4.9e-3), 2402 C(200,0),
2320 CMPLX(0,-5.1e-3), 2403 C(0,4.9e-3),
2321 CMPLX(0,2e-6), 2404 C(0,-5.1e-3),
2322 CMPLX(0,-2), 2405 C(0,2e-6),
2323 CMPLX(0,20), 2406 C(0,-2),
2324 CMPLX(0,-200), 2407 C(0,20),
2325 CMPLX(Inf,0), 2408 C(0,-200),
2326 CMPLX(-Inf,0), 2409 C(Inf,0),
2327 CMPLX(0,Inf), 2410 C(-Inf,0),
2328 CMPLX(0,-Inf), 2411 C(0,Inf),
2329 CMPLX(Inf,Inf), 2412 C(0,-Inf),
2330 CMPLX(Inf,-Inf), 2413 C(Inf,Inf),
2331 CMPLX(NaN,NaN), 2414 C(Inf,-Inf),
2332 CMPLX(NaN,0), 2415 C(NaN,NaN),
2333 CMPLX(0,NaN), 2416 C(NaN,0),
2334 CMPLX(NaN,Inf), 2417 C(0,NaN),
2335 CMPLX(Inf,NaN), 2418 C(NaN,Inf),
2336 CMPLX(39, 6.4e-5), 2419 C(Inf,NaN),
2337 CMPLX(41, 6.09e-5), 2420 C(39, 6.4e-5),
2338 CMPLX(4.9e7, 5e-11), 2421 C(41, 6.09e-5),
2339 CMPLX(5.1e7, 4.8e-11), 2422 C(4.9e7, 5e-11),
2340 CMPLX(1e9, 2.4e-12), 2423 C(5.1e7, 4.8e-11),
2341 CMPLX(1e11, 2.4e-14), 2424 C(1e9, 2.4e-12),
2342 CMPLX(1e13, 2.4e-16), 2425 C(1e11, 2.4e-14),
2343 CMPLX(1e300, 2.4e-303) 2426 C(1e13, 2.4e-16),
2427 C(1e300, 2.4e-303)
2344 }; 2428 };
2345 cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple 2429 cmplx w[NTST] = { // dawson(z[i]), evaluated with Maple
2346 CMPLX(0.1635394094345355614904345232875688576839, 2430 C(0.1635394094345355614904345232875688576839,
2347 -0.1531245755371229803585918112683241066853), 2431 -0.1531245755371229803585918112683241066853),
2348 CMPLX(-0.1635394094345355614904345232875688576839, 2432 C(-0.1635394094345355614904345232875688576839,
2349 -0.1531245755371229803585918112683241066853), 2433 -0.1531245755371229803585918112683241066853),
2350 CMPLX(0.1635394094345355614904345232875688576839, 2434 C(0.1635394094345355614904345232875688576839,
2351 0.1531245755371229803585918112683241066853), 2435 0.1531245755371229803585918112683241066853),
2352 CMPLX(-0.1635394094345355614904345232875688576839, 2436 C(-0.1635394094345355614904345232875688576839,
2353 0.1531245755371229803585918112683241066853), 2437 0.1531245755371229803585918112683241066853),
2354 CMPLX(-0.01619082256681596362895875232699626384420, 2438 C(-0.01619082256681596362895875232699626384420,
2355 -0.005210224203359059109181555401330902819419), 2439 -0.005210224203359059109181555401330902819419),
2356 CMPLX(0.01078377080978103125464543240346760257008, 2440 C(0.01078377080978103125464543240346760257008,
2357 0.006866888783433775382193630944275682670599), 2441 0.006866888783433775382193630944275682670599),
2358 CMPLX(-0.5808616819196736225612296471081337245459, 2442 C(-0.5808616819196736225612296471081337245459,
2359 0.6688593905505562263387760667171706325749), 2443 0.6688593905505562263387760667171706325749),
2360 CMPLX(Inf, 2444 C(Inf,
2361 -Inf), 2445 -Inf),
2362 CMPLX(0.1000052020902036118082966385855563526705e-7, 2446 C(0.1000052020902036118082966385855563526705e-7,
2363 0.005100088434920073153418834680320146441685), 2447 0.005100088434920073153418834680320146441685),
2364 CMPLX(0.004950156837581592745389973960217444687524, 2448 C(0.004950156837581592745389973960217444687524,
2365 -0.004899838305155226382584756154100963570500), 2449 -0.004899838305155226382584756154100963570500),
2366 CMPLX(0.4244534840871830045021143490355372016428, 2450 C(0.005100176864319675957314822982399286703798,
2367 0.002820278933186814021399602648373095266538), 2451 0.005099823128319785355949825238269336481254),
2368 CMPLX(-0.1021340733271046543881236523269967674156, 2452 C(0.4244534840871830045021143490355372016428,
2369 -0.00001045696456072005761498961861088944159916), 2453 0.002820278933186814021399602648373095266538),
2370 CMPLX(-0.01000200120119206748855061636187197886859, 2454 C(-0.1021340733271046543881236523269967674156,
2371 0.9805885888237419500266621041508714123763e-8), 2455 -0.00001045696456072005761498961861088944159916),
2372 CMPLX(0.001000002000012000023960527532953151819595, 2456 C(-0.01000200120119206748855061636187197886859,
2373 -0.9800058800588007290937355024646722133204e-11), 2457 0.9805885888237419500266621041508714123763e-8),
2374 CMPLX(0.4244549085628511778373438768121222815752, 2458 C(0.001000002000012000023960527532953151819595,
2375 0.002935393851311701428647152230552122898291), 2459 -0.9800058800588007290937355024646722133204e-11),
2376 CMPLX(-0.1021340732357117208743299813648493928105, 2460 C(0.4244549085628511778373438768121222815752,
2377 -0.00001088377943049851799938998805451564893540), 2461 0.002935393851311701428647152230552122898291),
2378 CMPLX(-0.01000200120119126652710792390331206563616, 2462 C(-0.1021340732357117208743299813648493928105,
2379 0.1020612612857282306892368985525393707486e-7), 2463 -0.00001088377943049851799938998805451564893540),
2380 CMPLX(0.1000000000007333333333344266666666664457e-5, 2464 C(-0.01000200120119126652710792390331206563616,
2381 0.2000000000001333333333323199999999978819e-5), 2465 0.1020612612857282306892368985525393707486e-7),
2382 CMPLX(0.1999999999994666666666675199999999990248e-5, 2466 C(0.1000000000007333333333344266666666664457e-5,
2383 0), 2467 0.2000000000001333333333323199999999978819e-5),
2384 CMPLX(0.3013403889237919660346644392864226952119, 2468 C(0.1999999999994666666666675199999999990248e-5,
2385 0), 2469 0),
2386 CMPLX(0.02503136792640367194699495234782353186858, 2470 C(0.3013403889237919660346644392864226952119,
2387 0), 2471 0),
2388 CMPLX(0.002500031251171948248596912483183760683918, 2472 C(0.02503136792640367194699495234782353186858,
2389 0), 2473 0),
2390 CMPLX(0,0.004900078433419939164774792850907128053308), 2474 C(0.002500031251171948248596912483183760683918,
2391 CMPLX(0,-0.005100088434920074173454208832365950009419), 2475 0),
2392 CMPLX(0,0.2000000000005333333333341866666666676419e-5), 2476 C(0,0.004900078433419939164774792850907128053308),
2393 CMPLX(0,-48.16001211429122974789822893525016528191), 2477 C(0,-0.005100088434920074173454208832365950009419),
2394 CMPLX(0,0.4627407029504443513654142715903005954668e174), 2478 C(0,0.2000000000005333333333341866666666676419e-5),
2395 CMPLX(0,-Inf), 2479 C(0,-48.16001211429122974789822893525016528191),
2396 CMPLX(0,0), 2480 C(0,0.4627407029504443513654142715903005954668e174),
2397 CMPLX(-0,0), 2481 C(0,-Inf),
2398 CMPLX(0, Inf), 2482 C(0,0),
2399 CMPLX(0, -Inf), 2483 C(-0,0),
2400 CMPLX(NaN, NaN), 2484 C(0, Inf),
2401 CMPLX(NaN, NaN), 2485 C(0, -Inf),
2402 CMPLX(NaN, NaN), 2486 C(NaN, NaN),
2403 CMPLX(NaN, 0), 2487 C(NaN, NaN),
2404 CMPLX(0, NaN), 2488 C(NaN, NaN),
2405 CMPLX(NaN, NaN), 2489 C(NaN, 0),
2406 CMPLX(NaN, NaN), 2490 C(0, NaN),
2407 CMPLX(0.01282473148489433743567240624939698290584, 2491 C(NaN, NaN),
2408 -0.2105957276516618621447832572909153498104e-7), 2492 C(NaN, NaN),
2409 CMPLX(0.01219875253423634378984109995893708152885, 2493 C(0.01282473148489433743567240624939698290584,
2410 -0.1813040560401824664088425926165834355953e-7), 2494 -0.2105957276516618621447832572909153498104e-7),
2411 CMPLX(0.1020408163265306334945473399689037886997e-7, 2495 C(0.01219875253423634378984109995893708152885,
2412 -0.1041232819658476285651490827866174985330e-25), 2496 -0.1813040560401824664088425926165834355953e-7),
2413 CMPLX(0.9803921568627452865036825956835185367356e-8, 2497 C(0.1020408163265306334945473399689037886997e-7,
2414 -0.9227220299884665067601095648451913375754e-26), 2498 -0.1041232819658476285651490827866174985330e-25),
2415 CMPLX(0.5000000000000000002500000000000000003750e-9, 2499 C(0.9803921568627452865036825956835185367356e-8,
2416 -0.1200000000000000001800000188712838420241e-29), 2500 -0.9227220299884665067601095648451913375754e-26),
2417 CMPLX(5.00000000000000000000025000000000000000000003e-12, 2501 C(0.5000000000000000002500000000000000003750e-9,
2418 -1.20000000000000000000018000000000000000000004e-36), 2502 -0.1200000000000000001800000188712838420241e-29),
2419 CMPLX(5.00000000000000000000000002500000000000000000e-14, 2503 C(5.00000000000000000000025000000000000000000003e-12,
2420 -1.20000000000000000000000001800000000000000000e-42), 2504 -1.20000000000000000000018000000000000000000004e-36),
2421 CMPLX(5e-301, 0) 2505 C(5.00000000000000000000000002500000000000000000e-14,
2506 -1.20000000000000000000000001800000000000000000e-42),
2507 C(5e-301, 0)
2422 }; 2508 };
2423 TST(Dawson, 1e-20); 2509 TST(Dawson, 1e-20);
2424 } 2510 }
2425 printf("#####################################\n"); 2511 printf("#####################################\n");
2426 printf("SUCCESS (max relative error = %g)\n", errmax_all); 2512 printf("SUCCESS (max relative error = %g)\n", errmax_all);