1/* Library libcerf: 2 * Compute complex error functions, based on a new implementation of 3 * Faddeeva's w_of_z. Also provide Dawson and Voigt functions. 4 * 5 * File err_fcts.c: 6 * Computate Dawson, Voigt, and several error functions, 7 * based on erfcx, im_w_of_x, w_of_z as implemented in separate files. 8 * 9 * Given w(z), the error functions are mostly straightforward 10 * to compute, except for certain regions where we have to 11 * switch to Taylor expansions to avoid cancellation errors 12 * [e.g. near the origin for erf(z)]. 13 * 14 * Copyright: 15 * (C) 2012 Massachusetts Institute of Technology 16 * (C) 2013 Forschungszentrum Jülich GmbH 17 * 18 * Licence: 19 * Permission is hereby granted, free of charge, to any person obtaining 20 * a copy of this software and associated documentation files (the 21 * "Software"), to deal in the Software without restriction, including 22 * without limitation the rights to use, copy, modify, merge, publish, 23 * distribute, sublicense, and/or sell copies of the Software, and to 24 * permit persons to whom the Software is furnished to do so, subject to 25 * the following conditions: 26 * 27 * The above copyright notice and this permission notice shall be 28 * included in all copies or substantial portions of the Software. 29 * 30 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 31 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 32 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 33 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 34 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 35 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 36 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 37 * 38 * Authors: 39 * Steven G. Johnson, Massachusetts Institute of Technology, 2012, core author 40 * Joachim Wuttke, Forschungszentrum Jülich, 2013, package maintainer 41 * 42 * Website: 43 * http://apps.jcns.fz-juelich.de/libcerf 44 * 45 * Revision history: 46 * ../CHANGELOG 47 * 48 * Man pages: 49 * cerf(3), dawson(3), voigt(3) 50 */ 51 52#include "cerf.h" 53#include <math.h> 54#include "defs.h" // defines _cerf_cmplx, NaN, C, cexp, ... 55 56const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 57const double s2pi = 2.5066282746310005024157652848110; // sqrt(2*pi) 58const double pi = 3.141592653589793238462643383279503; 59 60/******************************************************************************/ 61/* Simple wrappers: cerfcx, cerfi, erfi, dawson */ 62/******************************************************************************/ 63 64_cerf_cmplx cerfcx(_cerf_cmplx z) 65{ 66 // Compute erfcx(z) = exp(z^2) erfc(z), 67 // the complex underflow-compensated complementary error function, 68 // trivially related to Faddeeva's w_of_z. 69 70 return w_of_z(C(-cimag(z), creal(z))); 71} 72 73_cerf_cmplx cerfi(_cerf_cmplx z) 74{ 75 // Compute erfi(z) = -i erf(iz), 76 // the rotated complex error function. 77 78 _cerf_cmplx e = cerf(C(-cimag(z),creal(z))); 79 return C(cimag(e), -creal(e)); 80} 81 82double erfi(double x) 83{ 84 // Compute erfi(x) = -i erf(ix), 85 // the imaginary error function. 86 87 return x*x > 720 ? (x > 0 ? Inf : -Inf) : exp(x*x) * im_w_of_x(x); 88} 89 90double dawson(double x) 91{ 92 93 // Compute dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x), 94 // Dawson's integral for a real argument. 95 96 return spi2 * im_w_of_x(x); 97} 98 99double re_w_of_z( double x, double y ) 100{ 101 return creal( w_of_z( C(x,y) ) ); 102} 103 104double im_w_of_z( double x, double y ) 105{ 106 return cimag( w_of_z( C(x,y) ) ); 107} 108 109/******************************************************************************/ 110/* voigt */ 111/******************************************************************************/ 112 113double voigt( double x, double sigma, double gamma ) 114{ 115 // Joachim Wuttke, January 2013. 116 117 // Compute Voigt's convolution of a Gaussian 118 // G(x,sigma) = 1/sqrt(2*pi)/|sigma| * exp(-x^2/2/sigma^2) 119 // and a Lorentzian 120 // L(x,gamma) = |gamma| / pi / ( x^2 + gamma^2 ), 121 // namely 122 // voigt(x,sigma,gamma) = 123 // \int_{-infty}^{infty} dx' G(x',sigma) L(x-x',gamma) 124 // using the relation 125 // voigt(x,sigma,gamma) = Re{ w(z) } / sqrt(2*pi) / |sigma| 126 // with 127 // z = (x+i*|gamma|) / sqrt(2) / |sigma|. 128 129 // Reference: Abramowitz&Stegun (1964), formula (7.4.13). 130 131 double gam = gamma < 0 ? -gamma : gamma; 132 double sig = sigma < 0 ? -sigma : sigma; 133 134 if ( gam==0 ) { 135 if ( sig==0 ) { 136 // It's kind of a delta function 137 return x ? 0 : Inf; 138 } else { 139 // It's a pure Gaussian 140 return exp( -x*x/2/(sig*sig) ) / s2pi / sig; 141 } 142 } else { 143 if ( sig==0 ) { 144 // It's a pure Lorentzian 145 return gam / pi / (x*x + gam*gam); 146 } else { 147 // Regular case, both parameters are nonzero 148 _cerf_cmplx z = complex_mul_cr(C(x, gam), 1. / sqrt(2) / sig); 149 return creal( w_of_z(z) ) / s2pi / sig; 150 // TODO: correct and activate the following: 151// double w = sqrt(gam*gam+sig*sig); // to work in reduced units 152// _cerf_cmplx z = C(x/w,gam/w) / sqrt(2) / (sig/w); 153// return creal( w_of_z(z) ) / s2pi / (sig/w); 154 } 155 } 156} 157 158/******************************************************************************/ 159/* cerf */ 160/******************************************************************************/ 161 162_cerf_cmplx cerf(_cerf_cmplx z) 163{ 164 165 // Steven G. Johnson, October 2012. 166 167 // Compute erf(z), the complex error function, 168 // using w_of_z except for certain regions. 169 170 double x = creal(z), y = cimag(z); 171 172 if (y == 0) 173 return C(erf(x), y); // preserve sign of 0 174 if (x == 0) // handle separately for speed & handling of y = Inf or NaN 175 return C(x, // preserve sign of 0 176 /* handle y -> Inf limit manually, since 177 exp(y^2) -> Inf but Im[w(y)] -> 0, so 178 IEEE will give us a NaN when it should be Inf */ 179 y*y > 720 ? (y > 0 ? Inf : -Inf) 180 : exp(y*y) * im_w_of_x(y)); 181 182 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 183 double mIm_z2 = -2*x*y; // Im(-z^2) 184 if (mRe_z2 < -750) // underflow 185 return (x >= 0 ? C(1.0, 0.0) : C(-1.0, 0.0));; 186 187 /* Handle positive and negative x via different formulas, 188 using the mirror symmetries of w, to avoid overflow/underflow 189 problems from multiplying exponentially large and small quantities. */ 190 if (x >= 0) { 191 if (x < 8e-2) { 192 if (fabs(y) < 1e-2) 193 goto taylor; 194 else if (fabs(mIm_z2) < 5e-3 && x < 5e-3) 195 goto taylor_erfi; 196 } 197 /* don't use complex exp function, since that will produce spurious NaN 198 values when multiplying w in an overflow situation. */ 199 return complex_sub_rc(1.0, complex_mul_rc(exp(mRe_z2), complex_mul_cc(C(cos(mIm_z2), sin(mIm_z2)), w_of_z(C(-y, x))))); 200 } 201 else { // x < 0 202 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call 203 if (fabs(y) < 1e-2) 204 goto taylor; 205 else if (fabs(mIm_z2) < 5e-3 && x > -5e-3) 206 goto taylor_erfi; 207 } 208 else if (isnan(x)) 209 return C(NaN, y == 0 ? 0 : NaN); 210 /* don't use complex exp function, since that will produce spurious NaN 211 values when multiplying w in an overflow situation. */ 212 return complex_add_rc(-1.0, complex_mul_rc(exp(mRe_z2), complex_mul_cc(C(cos(mIm_z2), sin(mIm_z2)), w_of_z(C(y, -x))))); 213 214 } 215 216 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 217 // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...) 218taylor: 219 { 220 _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 221 return 222 complex_mul_cc(z, complex_add_rc(1.1283791670955125739, 223 complex_mul_cc(mz2, complex_add_rc(0.37612638903183752464, 224 complex_mul_cc(mz2, complex_add_rc(0.11283791670955125739, 225 complex_mul_cc(mz2, complex_add_rc(0.026866170645131251760, 226 complex_mul_cr(mz2, 0.0052239776254421878422))))))))); 227 228 229 } 230 231 /* for small |x| and small |xy|, 232 use Taylor series to avoid cancellation inaccuracy: 233 erf(x+iy) = erf(iy) 234 + 2*exp(y^2)/sqrt(pi) * 235 [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... 236 - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] 237 where: 238 erf(iy) = exp(y^2) * Im[w(y)] 239 */ 240taylor_erfi: 241 { 242 double x2 = x*x, y2 = y*y; 243 double expy2 = exp(y2); 244 return C 245 (expy2 * x * (1.1283791670955125739 246 - x2 * (0.37612638903183752464 247 + 0.75225277806367504925*y2) 248 + x2*x2 * (0.11283791670955125739 249 + y2 * (0.45135166683820502956 250 + 0.15045055561273500986*y2))), 251 expy2 * (im_w_of_x(y) 252 - x2*y * (1.1283791670955125739 253 - x2 * (0.56418958354775628695 254 + 0.37612638903183752464*y2)))); 255 } 256} // cerf 257 258/******************************************************************************/ 259/* cerfc */ 260/******************************************************************************/ 261 262_cerf_cmplx cerfc(_cerf_cmplx z) 263{ 264 // Steven G. Johnson, October 2012. 265 266 // Compute erfc(z) = 1 - erf(z), the complex complementary error function, 267 // using w_of_z except for certain regions. 268 269 double x = creal(z), y = cimag(z); 270 271 if (x == 0.) 272 return C(1, 273 /* handle y -> Inf limit manually, since 274 exp(y^2) -> Inf but Im[w(y)] -> 0, so 275 IEEE will give us a NaN when it should be Inf */ 276 y*y > 720 ? (y > 0 ? -Inf : Inf) 277 : -exp(y*y) * im_w_of_x(y)); 278 if (y == 0.) { 279 if (x*x > 750) // underflow 280 return C(x >= 0 ? 0.0 : 2.0, 281 -y); // preserve sign of 0 282 return C(x >= 0 ? exp(-x*x) * erfcx(x) 283 : 2. - exp(-x*x) * erfcx(-x), 284 -y); // preserve sign of zero 285 } 286 287 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 288 double mIm_z2 = -2*x*y; // Im(-z^2) 289 if (mRe_z2 < -750) // underflow 290 return C((x >= 0 ? 0.0 : 2.0), 0.0); 291 292 if (x >= 0) 293 return cexp(complex_mul_cc(C(mRe_z2, mIm_z2), w_of_z(C(-y,x)))); 294 else 295 return complex_sub_rc(2.0, complex_mul_cc(cexp(C(mRe_z2, mIm_z2)), w_of_z(C(y, -x)))); 296} // cerfc 297 298/******************************************************************************/ 299/* cdawson */ 300/******************************************************************************/ 301 302_cerf_cmplx cdawson(_cerf_cmplx z) 303{ 304 305 // Steven G. Johnson, October 2012. 306 307 // Compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z), 308 // Dawson's integral for a complex argument, 309 // using w_of_z except for certain regions. 310 311 double x = creal(z), y = cimag(z); 312 313 // handle axes separately for speed & proper handling of x or y = Inf or NaN 314 if (y == 0) 315 return C(spi2 * im_w_of_x(x), 316 -y); // preserve sign of 0 317 if (x == 0) { 318 double y2 = y*y; 319 if (y2 < 2.5e-5) { // Taylor expansion 320 return C(x, // preserve sign of 0 321 y * (1. 322 + y2 * (0.6666666666666666666666666666666666666667 323 + y2 * 0.26666666666666666666666666666666666667))); 324 } 325 return C(x, // preserve sign of 0 326 spi2 * (y >= 0 327 ? exp(y2) - erfcx(y) 328 : erfcx(-y) - exp(y2))); 329 } 330 331 double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow 332 double mIm_z2 = -2*x*y; // Im(-z^2) 333 _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 334 335 /* Handle positive and negative x via different formulas, 336 using the mirror symmetries of w, to avoid overflow/underflow 337 problems from multiplying exponentially large and small quantities. */ 338 if (y >= 0) { 339 if (y < 5e-3) { 340 if (fabs(x) < 5e-3) 341 goto taylor; 342 else if (fabs(mIm_z2) < 5e-3) 343 goto taylor_realaxis; 344 } 345 _cerf_cmplx res = complex_sub_cc(cexp(mz2), w_of_z(z)); 346 return complex_mul_rc(spi2, C(-cimag(res), creal(res))); 347 } 348 else { // y < 0 349 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call 350 if (fabs(x) < 5e-3) 351 goto taylor; 352 else if (fabs(mIm_z2) < 5e-3) 353 goto taylor_realaxis; 354 } 355 else if (isnan(y)) 356 return C(x == 0 ? 0 : NaN, NaN); 357 { 358 _cerf_cmplx res = complex_sub_cc(w_of_z(complex_neg(z)), cexp(mz2)); 359 return complex_mul_rc(spi2, C(-cimag(res), creal(res))); 360 } 361 } 362 363 // Use Taylor series for small |z|, to avoid cancellation inaccuracy 364 // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... 365taylor: 366 return complex_mul_cc(z, complex_add_rc(1., 367 complex_mul_cc(mz2, complex_add_rc(0.6666666666666666666666666666666666666667, 368 complex_mul_cr(mz2, 0.2666666666666666666666666666666666666667))))); 369 /* for small |y| and small |xy|, 370 use Taylor series to avoid cancellation inaccuracy: 371 dawson(x + iy) 372 = D + y^2 (D + x - 2Dx^2) 373 + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3) 374 + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3) 375 + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ... 376 where D = dawson(x) 377 378 However, for large |x|, 2Dx -> 1 which gives cancellation problems in 379 this series (many of the leading terms cancel). So, for large |x|, 380 we need to substitute a continued-fraction expansion for D. 381 382 dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...)))))) 383 384 The 6 terms shown here seems to be the minimum needed to be 385 accurate as soon as the simpler Taylor expansion above starts 386 breaking down. Using this 6-term expansion, factoring out the 387 denominator, and simplifying with Maple, we obtain: 388 389 Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x 390 = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4 391 Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y 392 = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4 393 394 Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction 395 expansion for the real part, and a 2-term expansion for the imaginary 396 part. (This avoids overflow problems for huge |x|.) This yields: 397 398 Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x) 399 Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1) 400 401 */ 402taylor_realaxis: 403 { 404 double x2 = x*x; 405 if (x2 > 1600) { // |x| > 40 406 double y2 = y*y; 407 if (x2 > 25e14) {// |x| > 5e7 408 double xy2 = (x*y)*(x*y); 409 return C((0.5 + y2 * (0.5 + 0.25*y2 410 - 0.16666666666666666667*xy2)) / x, 411 y * (-1 + y2 * (-0.66666666666666666667 412 + 0.13333333333333333333*xy2 413 - 0.26666666666666666667*y2)) 414 / (2*x2 - 1)); 415 } 416 return complex_mul_rc((1. / (-15 + x2 * (90 + x2 * (-60 + 8 * x2)))), 417 C(x * (33 + x2 * (-28 + 4 * x2) 418 + +y2 * (18 - 4 * x2 + 4 * y2)), 419 +y * (-15 + x2 * (24 - 4 * x2) 420 + +y2 * (4 * x2 - 10 - 4 * y2)))); 421 } 422 else { 423 double D = spi2 * im_w_of_x(x); 424 double y2 = y*y; 425 return C 426 (D + y2 * (D + x - 2*D*x2) 427 + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2)) 428 + x * (0.83333333333333333333 429 - 0.33333333333333333333 * x2)), 430 y * (1 - 2*D*x 431 + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2)) 432 + y2*y2 * (0.26666666666666666667 - 433 x2 * (0.6 - 0.13333333333333333333 * x2) 434 - D*x * (1 - x2 * (1.3333333333333333333 435 - 0.26666666666666666667 * x2))))); 436 } 437 } 438} // cdawson 439 |