/* Library libcerf: * Compute complex error functions, based on a new implementation of * Faddeeva's w_of_z. Also provide Dawson and Voigt functions. * * File w_of_z.c: * Computation of Faddeeva's complex scaled error function, * w(z) = exp(-z^2) * erfc(-i*z), * nameless function (7.1.3) of Abramowitz&Stegun (1964), * also known as the plasma dispersion function. * * This implementation uses a combination of different algorithms. * See man 3 w_of_z for references. * * 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 page: * w_of_z(3) */ /* Todo: use local declarations (older compilers) (HH). */ /* Computes various error functions (erf, erfc, erfi, erfcx), including the Dawson integral, in the complex plane, based on algorithms for the computation of the Faddeeva function w(z) = exp(-z^2) * erfc(-i*z). 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)]. */ #include "cerf.h" #include #include #include "defs.h" // defines _cerf_cmplx, NaN, C, cexp, ... // for analysing the algorithm: EXPORT int faddeeva_algorithm; EXPORT int faddeeva_nofterms; /******************************************************************************/ /* auxiliary functions */ /******************************************************************************/ static inline double sinc(double x, double sinx) { // return sinc(x) = sin(x)/x, given both x and sin(x) // [since we only use this in cases where sin(x) has already been computed] return fabs(x) < 1e-4 ? 1 - (0.1666666666666666666667)*x*x : sinx / x; } static inline double sinh_taylor(double x) { // sinh(x) via Taylor series, accurate to machine precision for |x| < 1e-2 return x * (1 + (x*x) * (0.1666666666666666666667 + 0.00833333333333333333333 * (x*x))); } static inline double sqr(double x) { return x*x; } /******************************************************************************/ /* precomputed table of expa2n2[n-1] = exp(-a2*n*n) */ /* for double-precision a2 = 0.26865... in w_of_z, below. */ /******************************************************************************/ static const double expa2n2[] = { 7.64405281671221563e-01, 3.41424527166548425e-01, 8.91072646929412548e-02, 1.35887299055460086e-02, 1.21085455253437481e-03, 6.30452613933449404e-05, 1.91805156577114683e-06, 3.40969447714832381e-08, 3.54175089099469393e-10, 2.14965079583260682e-12, 7.62368911833724354e-15, 1.57982797110681093e-17, 1.91294189103582677e-20, 1.35344656764205340e-23, 5.59535712428588720e-27, 1.35164257972401769e-30, 1.90784582843501167e-34, 1.57351920291442930e-38, 7.58312432328032845e-43, 2.13536275438697082e-47, 3.51352063787195769e-52, 3.37800830266396920e-57, 1.89769439468301000e-62, 6.22929926072668851e-68, 1.19481172006938722e-73, 1.33908181133005953e-79, 8.76924303483223939e-86, 3.35555576166254986e-92, 7.50264110688173024e-99, 9.80192200745410268e-106, 7.48265412822268959e-113, 3.33770122566809425e-120, 8.69934598159861140e-128, 1.32486951484088852e-135, 1.17898144201315253e-143, 6.13039120236180012e-152, 1.86258785950822098e-160, 3.30668408201432783e-169, 3.43017280887946235e-178, 2.07915397775808219e-187, 7.36384545323984966e-197, 1.52394760394085741e-206, 1.84281935046532100e-216, 1.30209553802992923e-226, 5.37588903521080531e-237, 1.29689584599763145e-247, 1.82813078022866562e-258, 1.50576355348684241e-269, 7.24692320799294194e-281, 2.03797051314726829e-292, 3.34880215927873807e-304, 0.0 // underflow (also prevents reads past array end, below) }; // expa2n2 /******************************************************************************/ /* w_of_z, Faddeeva's scaled complex error function */ /******************************************************************************/ _cerf_cmplx w_of_z(_cerf_cmplx z) { faddeeva_nofterms = 0; // Steven G. Johnson, October 2012. if (creal(z) == 0.0) { // Purely imaginary input, purely real output. // However, use creal(z) to give correct sign of 0 in cimag(w). return C(erfcx(cimag(z)), creal(z)); } if (cimag(z) == 0) { // Purely real input, complex output. return C(exp(-sqr(creal(z))), im_w_of_x(creal(z))); } const double relerr = DBL_EPSILON; const double a = 0.518321480430085929872; // pi / sqrt(-log(eps*0.5)) const double c = 0.329973702884629072537; // (2/pi) * a; const double a2 = 0.268657157075235951582; // a^2 const double x = fabs(creal(z)); const double y = cimag(z); const double ya = fabs(y); _cerf_cmplx ret = C(0., 0.); // return value double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0, sum5 = 0; if (ya > 7 || (x > 6 // continued fraction is faster /* As pointed out by M. Zaghloul, the continued fraction seems to give a large relative error in Re w(z) for |x| ~ 6 and small |y|, so use algorithm 816 in this region: */ && (ya > 0.1 || (x > 8 && ya > 1e-10) || x > 28))) { faddeeva_algorithm = 100; /* Poppe & Wijers suggest using a number of terms nu = 3 + 1442 / (26*rho + 77) where rho = sqrt((x/x0)^2 + (y/y0)^2) where x0=6.3, y0=4.4. (They only use this expansion for rho >= 1, but rho a little less than 1 seems okay too.) Instead, I did my own fit to a slightly different function that avoids the hypotenuse calculation, using NLopt to minimize the sum of the squares of the errors in nu with the constraint that the estimated nu be >= minimum nu to attain machine precision. I also separate the regions where nu == 2 and nu == 1. */ const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi) double xs = y < 0 ? -creal(z) : creal(z); // compute for -z if y < 0 if (x + ya > 4000) { // nu <= 2 if (x + ya > 1e7) { // nu == 1, w(z) = i/sqrt(pi) / z // scale to avoid overflow if (x > ya) { faddeeva_algorithm += 1; double yax = ya / xs; faddeeva_algorithm = 100; double denom = ispi / (xs + yax*ya); ret = C(denom*yax, denom); } else if (isinf(ya)) { faddeeva_algorithm += 2; return ((isnan(x) || y < 0) ? C(NaN,NaN) : C(0,0)); } else { faddeeva_algorithm += 3; double xya = xs / ya; double denom = ispi / (xya*xs + ya); ret = C(denom, denom*xya); } } else { // nu == 2, w(z) = i/sqrt(pi) * z / (z*z - 0.5) faddeeva_algorithm += 4; double dr = xs*xs - ya*ya - 0.5, di = 2*xs*ya; double denom = ispi / (dr*dr + di*di); ret = C(denom * (xs*di-ya*dr), denom * (xs*dr+ya*di)); } } else { // compute nu(z) estimate and do general continued fraction faddeeva_algorithm += 5; const double c0=3.9, c1=11.398, c2=0.08254, c3=0.1421, c4=0.2023; // fit double nu = floor(c0 + c1 / (c2*x + c3*ya + c4)); double wr = xs, wi = ya; for (nu = 0.5 * (nu - 1); nu > 0.4; nu -= 0.5) { // w <- z - nu/w: double denom = nu / (wr*wr + wi*wi); wr = xs - wr * denom; wi = ya + wi * denom; } { // w(z) = i/sqrt(pi) / w: double denom = ispi / (wr*wr + wi*wi); ret = C(denom*wi, denom*wr); } } if (y < 0) { faddeeva_algorithm += 10; // use w(z) = 2.0*exp(-z*z) - w(-z), // but be careful of overflow in exp(-z*z) // = exp(-(xs*xs-ya*ya) -2*i*xs*ya) return complex_sub_cc(complex_mul_rc(2.0,cexp(C((ya-xs)*(xs+ya), 2*xs*y))), ret); } else return ret; } /* Note: The test that seems to be suggested in the paper is x < sqrt(-log(DBL_MIN)), about 26.6, since otherwise exp(-x^2) underflows to zero and sum1,sum2,sum4 are zero. However, long before this occurs, the sum1,sum2,sum4 contributions are negligible in double precision; I find that this happens for x > about 6, for all y. On the other hand, I find that the case where we compute all of the sums is faster (at least with the precomputed expa2n2 table) until about x=10. Furthermore, if we try to compute all of the sums for x > 20, I find that we sometimes run into numerical problems because underflow/overflow problems start to appear in the various coefficients of the sums, below. Therefore, we use x < 10 here. */ else if (x < 10) { faddeeva_algorithm = 200; double prod2ax = 1, prodm2ax = 1; double expx2; if (isnan(y)) { faddeeva_algorithm += 99; return C(y,y); } if (x < 5e-4) { // compute sum4 and sum5 together as sum5-sum4 // This special case is needed for accuracy. faddeeva_algorithm += 1; const double x2 = x*x; expx2 = 1 - x2 * (1 - 0.5*x2); // exp(-x*x) via Taylor // compute exp(2*a*x) and exp(-2*a*x) via Taylor, to double precision const double ax2 = 1.036642960860171859744*x; // 2*a*x const double exp2ax = 1 + ax2 * (1 + ax2 * (0.5 + 0.166666666666666666667*ax2)); const double expm2ax = 1 - ax2 * (1 - ax2 * (0.5 - 0.166666666666666666667*ax2)); for (int n = 1; ; ++n) { ++faddeeva_nofterms; const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); prod2ax *= exp2ax; prodm2ax *= expm2ax; sum1 += coef; sum2 += coef * prodm2ax; sum3 += coef * prod2ax; // really = sum5 - sum4 sum5 += coef * (2*a) * n * sinh_taylor((2*a)*n*x); // test convergence via sum3 if (coef * prod2ax < relerr * sum3) break; } } else { // x > 5e-4, compute sum4 and sum5 separately faddeeva_algorithm += 2; expx2 = exp(-x*x); const double exp2ax = exp((2*a)*x), expm2ax = 1 / exp2ax; for (int n = 1; ; ++n) { ++faddeeva_nofterms; const double coef = expa2n2[n-1] * expx2 / (a2*(n*n) + y*y); prod2ax *= exp2ax; prodm2ax *= expm2ax; sum1 += coef; sum2 += coef * prodm2ax; sum4 += (coef * prodm2ax) * (a*n); sum3 += coef * prod2ax; sum5 += (coef * prod2ax) * (a*n); // test convergence via sum5, since this sum has the slowest decay if ((coef * prod2ax) * (a*n) < relerr * sum5) break; } } const double expx2erfcxy = // avoid spurious overflow for large negative y y > -6 // for y < -6, erfcx(y) = 2*exp(y*y) to double precision ? expx2*erfcx(y) : 2*exp(y*y-x*x); if (y > 5) { // imaginary terms cancel faddeeva_algorithm += 10; const double sinxy = sin(x*y); ret = C((expx2erfcxy - c*y*sum1) * cos(2*x*y) + (c*x*expx2) * sinxy * sinc(x*y, sinxy), 0.0); } else { faddeeva_algorithm += 20; double xs = creal(z); const double sinxy = sin(xs*y); const double sin2xy = sin(2*xs*y), cos2xy = cos(2*xs*y); const double coef1 = expx2erfcxy - c*y*sum1; const double coef2 = c*xs*expx2; ret = C(coef1 * cos2xy + coef2 * sinxy * sinc(xs*y, sinxy), coef2 * sinc(2*xs*y, sin2xy) - coef1 * sin2xy); } } else { // x large: only sum3 & sum5 contribute (see above note) faddeeva_algorithm = 300; if (isnan(x)) return C(x,x); if (isnan(y)) return C(y,y); ret = C(exp(-x*x),0.0); // |y| < 1e-10, so we only need exp(-x*x) term // (round instead of ceil as in original paper; note that x/a > 1 here) double n0 = floor(x/a + 0.5); // sum in both directions, starting at n0 double dx = a*n0 - x; sum3 = exp(-dx*dx) / (a2*(n0*n0) + y*y); sum5 = a*n0 * sum3; double exp1 = exp(4*a*dx), exp1dn = 1; int dn; for (dn = 1; n0 - dn > 0; ++dn) { // loop over n0-dn and n0+dn terms double np = n0 + dn, nm = n0 - dn; double tp = exp(-sqr(a*dn+dx)); double tm = tp * (exp1dn *= exp1); // trick to get tm from tp tp /= (a2*(np*np) + y*y); tm /= (a2*(nm*nm) + y*y); sum3 += tp + tm; sum5 += a * (np * tp + nm * tm); if (a * (np * tp + nm * tm) < relerr * sum5) goto finish; } while (1) { // loop over n0+dn terms only (since n0-dn <= 0) double np = n0 + dn++; double tp = exp(-sqr(a*dn+dx)) / (a2*(np*np) + y*y); sum3 += tp; sum5 += a * np * tp; if (a * np * tp < relerr * sum5) goto finish; } } finish: return complex_add_cc(ret, C((0.5*c)*y*(sum2+sum3),(0.5*c)*copysign(sum5-sum4, creal(z)))); } // w_of_z