/* Library libcerf: * Compute complex error functions, based on a new implementation of * Faddeeva's w_of_z. Also provide Dawson and Voigt functions. * * File err_fcts.c: * Computate Dawson, Voigt, and several error functions, * based on erfcx, im_w_of_x, w_of_z as implemented in separate files. * * Given w(z), the error functions are mostly straightforward * to compute, except for certain regions where we have to * switch to Taylor expansions to avoid cancellation errors * [e.g. near the origin for erf(z)]. * * Copyright: * (C) 2012 Massachusetts Institute of Technology * (C) 2013 Forschungszentrum Jülich GmbH * * Licence: * Permission is hereby granted, free of charge, to any person obtaining * a copy of this software and associated documentation files (the * "Software"), to deal in the Software without restriction, including * without limitation the rights to use, copy, modify, merge, publish, * distribute, sublicense, and/or sell copies of the Software, and to * permit persons to whom the Software is furnished to do so, subject to * the following conditions: * * The above copyright notice and this permission notice shall be * included in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * * Authors: * Steven G. Johnson, Massachusetts Institute of Technology, 2012, core author * Joachim Wuttke, Forschungszentrum Jülich, 2013, package maintainer * * Website: * http://apps.jcns.fz-juelich.de/libcerf * * Revision history: * ../CHANGELOG * * Man pages: * cerf(3), dawson(3), voigt(3) */ #include "cerf.h" #include #include "defs.h" // defines _cerf_cmplx, NaN, C, cexp, ... const double spi2 = 0.8862269254527580136490837416705725913990; // sqrt(pi)/2 const double s2pi = 2.5066282746310005024157652848110; // sqrt(2*pi) const double pi = 3.141592653589793238462643383279503; /******************************************************************************/ /* Simple wrappers: cerfcx, cerfi, erfi, dawson */ /******************************************************************************/ _cerf_cmplx cerfcx(_cerf_cmplx z) { // Compute erfcx(z) = exp(z^2) erfc(z), // the complex underflow-compensated complementary error function, // trivially related to Faddeeva's w_of_z. return w_of_z(C(-cimag(z), creal(z))); } _cerf_cmplx cerfi(_cerf_cmplx z) { // Compute erfi(z) = -i erf(iz), // the rotated complex error function. _cerf_cmplx e = cerf(C(-cimag(z),creal(z))); return C(cimag(e), -creal(e)); } double erfi(double x) { // Compute erfi(x) = -i erf(ix), // the imaginary error function. return x*x > 720 ? (x > 0 ? Inf : -Inf) : exp(x*x) * im_w_of_x(x); } double dawson(double x) { // Compute dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x), // Dawson's integral for a real argument. return spi2 * im_w_of_x(x); } double re_w_of_z( double x, double y ) { return creal( w_of_z( C(x,y) ) ); } double im_w_of_z( double x, double y ) { return cimag( w_of_z( C(x,y) ) ); } /******************************************************************************/ /* voigt */ /******************************************************************************/ double voigt( double x, double sigma, double gamma ) { // Joachim Wuttke, January 2013. // Compute Voigt's convolution of a Gaussian // G(x,sigma) = 1/sqrt(2*pi)/|sigma| * exp(-x^2/2/sigma^2) // and a Lorentzian // L(x,gamma) = |gamma| / pi / ( x^2 + gamma^2 ), // namely // voigt(x,sigma,gamma) = // \int_{-infty}^{infty} dx' G(x',sigma) L(x-x',gamma) // using the relation // voigt(x,sigma,gamma) = Re{ w(z) } / sqrt(2*pi) / |sigma| // with // z = (x+i*|gamma|) / sqrt(2) / |sigma|. // Reference: Abramowitz&Stegun (1964), formula (7.4.13). double gam = gamma < 0 ? -gamma : gamma; double sig = sigma < 0 ? -sigma : sigma; if ( gam==0 ) { if ( sig==0 ) { // It's kind of a delta function return x ? 0 : Inf; } else { // It's a pure Gaussian return exp( -x*x/2/(sig*sig) ) / s2pi / sig; } } else { if ( sig==0 ) { // It's a pure Lorentzian return gam / pi / (x*x + gam*gam); } else { // Regular case, both parameters are nonzero _cerf_cmplx z = complex_mul_cr(C(x, gam), 1. / sqrt(2) / sig); return creal( w_of_z(z) ) / s2pi / sig; // TODO: correct and activate the following: // double w = sqrt(gam*gam+sig*sig); // to work in reduced units // _cerf_cmplx z = C(x/w,gam/w) / sqrt(2) / (sig/w); // return creal( w_of_z(z) ) / s2pi / (sig/w); } } } /******************************************************************************/ /* cerf */ /******************************************************************************/ _cerf_cmplx cerf(_cerf_cmplx z) { // Steven G. Johnson, October 2012. // Compute erf(z), the complex error function, // using w_of_z except for certain regions. double x = creal(z), y = cimag(z); if (y == 0) return C(erf(x), y); // preserve sign of 0 if (x == 0) // handle separately for speed & handling of y = Inf or NaN return C(x, // preserve sign of 0 /* handle y -> Inf limit manually, since exp(y^2) -> Inf but Im[w(y)] -> 0, so IEEE will give us a NaN when it should be Inf */ y*y > 720 ? (y > 0 ? Inf : -Inf) : exp(y*y) * im_w_of_x(y)); double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow double mIm_z2 = -2*x*y; // Im(-z^2) if (mRe_z2 < -750) // underflow return (x >= 0 ? C(1.0, 0.0) : C(-1.0, 0.0));; /* Handle positive and negative x via different formulas, using the mirror symmetries of w, to avoid overflow/underflow problems from multiplying exponentially large and small quantities. */ if (x >= 0) { if (x < 8e-2) { if (fabs(y) < 1e-2) goto taylor; else if (fabs(mIm_z2) < 5e-3 && x < 5e-3) goto taylor_erfi; } /* don't use complex exp function, since that will produce spurious NaN values when multiplying w in an overflow situation. */ 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))))); } else { // x < 0 if (x > -8e-2) { // duplicate from above to avoid fabs(x) call if (fabs(y) < 1e-2) goto taylor; else if (fabs(mIm_z2) < 5e-3 && x > -5e-3) goto taylor_erfi; } else if (isnan(x)) return C(NaN, y == 0 ? 0 : NaN); /* don't use complex exp function, since that will produce spurious NaN values when multiplying w in an overflow situation. */ 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))))); } // Use Taylor series for small |z|, to avoid cancellation inaccuracy // erf(z) = 2/sqrt(pi) * z * (1 - z^2/3 + z^4/10 - z^6/42 + z^8/216 + ...) taylor: { _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 return complex_mul_cc(z, complex_add_rc(1.1283791670955125739, complex_mul_cc(mz2, complex_add_rc(0.37612638903183752464, complex_mul_cc(mz2, complex_add_rc(0.11283791670955125739, complex_mul_cc(mz2, complex_add_rc(0.026866170645131251760, complex_mul_cr(mz2, 0.0052239776254421878422))))))))); } /* for small |x| and small |xy|, use Taylor series to avoid cancellation inaccuracy: erf(x+iy) = erf(iy) + 2*exp(y^2)/sqrt(pi) * [ x * (1 - x^2 * (1+2y^2)/3 + x^4 * (3+12y^2+4y^4)/30 + ... - i * x^2 * y * (1 - x^2 * (3+2y^2)/6 + ...) ] where: erf(iy) = exp(y^2) * Im[w(y)] */ taylor_erfi: { double x2 = x*x, y2 = y*y; double expy2 = exp(y2); return C (expy2 * x * (1.1283791670955125739 - x2 * (0.37612638903183752464 + 0.75225277806367504925*y2) + x2*x2 * (0.11283791670955125739 + y2 * (0.45135166683820502956 + 0.15045055561273500986*y2))), expy2 * (im_w_of_x(y) - x2*y * (1.1283791670955125739 - x2 * (0.56418958354775628695 + 0.37612638903183752464*y2)))); } } // cerf /******************************************************************************/ /* cerfc */ /******************************************************************************/ _cerf_cmplx cerfc(_cerf_cmplx z) { // Steven G. Johnson, October 2012. // Compute erfc(z) = 1 - erf(z), the complex complementary error function, // using w_of_z except for certain regions. double x = creal(z), y = cimag(z); if (x == 0.) return C(1, /* handle y -> Inf limit manually, since exp(y^2) -> Inf but Im[w(y)] -> 0, so IEEE will give us a NaN when it should be Inf */ y*y > 720 ? (y > 0 ? -Inf : Inf) : -exp(y*y) * im_w_of_x(y)); if (y == 0.) { if (x*x > 750) // underflow return C(x >= 0 ? 0.0 : 2.0, -y); // preserve sign of 0 return C(x >= 0 ? exp(-x*x) * erfcx(x) : 2. - exp(-x*x) * erfcx(-x), -y); // preserve sign of zero } double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow double mIm_z2 = -2*x*y; // Im(-z^2) if (mRe_z2 < -750) // underflow return C((x >= 0 ? 0.0 : 2.0), 0.0); if (x >= 0) return cexp(complex_mul_cc(C(mRe_z2, mIm_z2), w_of_z(C(-y,x)))); else return complex_sub_rc(2.0, complex_mul_cc(cexp(C(mRe_z2, mIm_z2)), w_of_z(C(y, -x)))); } // cerfc /******************************************************************************/ /* cdawson */ /******************************************************************************/ _cerf_cmplx cdawson(_cerf_cmplx z) { // Steven G. Johnson, October 2012. // Compute Dawson(z) = sqrt(pi)/2 * exp(-z^2) * erfi(z), // Dawson's integral for a complex argument, // using w_of_z except for certain regions. double x = creal(z), y = cimag(z); // handle axes separately for speed & proper handling of x or y = Inf or NaN if (y == 0) return C(spi2 * im_w_of_x(x), -y); // preserve sign of 0 if (x == 0) { double y2 = y*y; if (y2 < 2.5e-5) { // Taylor expansion return C(x, // preserve sign of 0 y * (1. + y2 * (0.6666666666666666666666666666666666666667 + y2 * 0.26666666666666666666666666666666666667))); } return C(x, // preserve sign of 0 spi2 * (y >= 0 ? exp(y2) - erfcx(y) : erfcx(-y) - exp(y2))); } double mRe_z2 = (y - x) * (x + y); // Re(-z^2), being careful of overflow double mIm_z2 = -2*x*y; // Im(-z^2) _cerf_cmplx mz2 = C(mRe_z2, mIm_z2); // -z^2 /* Handle positive and negative x via different formulas, using the mirror symmetries of w, to avoid overflow/underflow problems from multiplying exponentially large and small quantities. */ if (y >= 0) { if (y < 5e-3) { if (fabs(x) < 5e-3) goto taylor; else if (fabs(mIm_z2) < 5e-3) goto taylor_realaxis; } _cerf_cmplx res = complex_sub_cc(cexp(mz2), w_of_z(z)); return complex_mul_rc(spi2, C(-cimag(res), creal(res))); } else { // y < 0 if (y > -5e-3) { // duplicate from above to avoid fabs(x) call if (fabs(x) < 5e-3) goto taylor; else if (fabs(mIm_z2) < 5e-3) goto taylor_realaxis; } else if (isnan(y)) return C(x == 0 ? 0 : NaN, NaN); { _cerf_cmplx res = complex_sub_cc(w_of_z(complex_neg(z)), cexp(mz2)); return complex_mul_rc(spi2, C(-cimag(res), creal(res))); } } // Use Taylor series for small |z|, to avoid cancellation inaccuracy // dawson(z) = z - 2/3 z^3 + 4/15 z^5 + ... taylor: return complex_mul_cc(z, complex_add_rc(1., complex_mul_cc(mz2, complex_add_rc(0.6666666666666666666666666666666666666667, complex_mul_cr(mz2, 0.2666666666666666666666666666666666666667))))); /* for small |y| and small |xy|, use Taylor series to avoid cancellation inaccuracy: dawson(x + iy) = D + y^2 (D + x - 2Dx^2) + y^4 (D/2 + 5x/6 - 2Dx^2 - x^3/3 + 2Dx^4/3) + iy [ (1-2Dx) + 2/3 y^2 (1 - 3Dx - x^2 + 2Dx^3) + y^4/15 (4 - 15Dx - 9x^2 + 20Dx^3 + 2x^4 - 4Dx^5) ] + ... where D = dawson(x) However, for large |x|, 2Dx -> 1 which gives cancellation problems in this series (many of the leading terms cancel). So, for large |x|, we need to substitute a continued-fraction expansion for D. dawson(x) = 0.5 / (x-0.5/(x-1/(x-1.5/(x-2/(x-2.5/(x...)))))) The 6 terms shown here seems to be the minimum needed to be accurate as soon as the simpler Taylor expansion above starts breaking down. Using this 6-term expansion, factoring out the denominator, and simplifying with Maple, we obtain: Re dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / x = 33 - 28x^2 + 4x^4 + y^2 (18 - 4x^2) + 4 y^4 Im dawson(x + iy) * (-15 + 90x^2 - 60x^4 + 8x^6) / y = -15 + 24x^2 - 4x^4 + 2/3 y^2 (6x^2 - 15) - 4 y^4 Finally, for |x| > 5e7, we can use a simpler 1-term continued-fraction expansion for the real part, and a 2-term expansion for the imaginary part. (This avoids overflow problems for huge |x|.) This yields: Re dawson(x + iy) = [1 + y^2 (1 + y^2/2 - (xy)^2/3)] / (2x) Im dawson(x + iy) = y [ -1 - 2/3 y^2 + y^4/15 (2x^2 - 4) ] / (2x^2 - 1) */ taylor_realaxis: { double x2 = x*x; if (x2 > 1600) { // |x| > 40 double y2 = y*y; if (x2 > 25e14) {// |x| > 5e7 double xy2 = (x*y)*(x*y); return C((0.5 + y2 * (0.5 + 0.25*y2 - 0.16666666666666666667*xy2)) / x, y * (-1 + y2 * (-0.66666666666666666667 + 0.13333333333333333333*xy2 - 0.26666666666666666667*y2)) / (2*x2 - 1)); } return complex_mul_rc((1. / (-15 + x2 * (90 + x2 * (-60 + 8 * x2)))), C(x * (33 + x2 * (-28 + 4 * x2) + +y2 * (18 - 4 * x2 + 4 * y2)), +y * (-15 + x2 * (24 - 4 * x2) + +y2 * (4 * x2 - 10 - 4 * y2)))); } else { double D = spi2 * im_w_of_x(x); double y2 = y*y; return C (D + y2 * (D + x - 2*D*x2) + y2*y2 * (D * (0.5 - x2 * (2 - 0.66666666666666666667*x2)) + x * (0.83333333333333333333 - 0.33333333333333333333 * x2)), y * (1 - 2*D*x + y2 * 0.66666666666666666667 * (1 - x2 - D*x * (3 - 2*x2)) + y2*y2 * (0.26666666666666666667 - x2 * (0.6 - 0.13333333333333333333 * x2) - D*x * (1 - x2 * (1.3333333333333333333 - 0.26666666666666666667 * x2))))); } } } // cdawson