w_of_z, im_w_of_x - Faddeeva's rescaled complex error function
#include <cerf.h>
double _Complex w_of_z ( double _Complex z );
double im_w_of_x ( double x );
Faddeeva's rescaled complex error function w(z), also called the plasma dispersion function.
w_of_z returns w(z) = exp(-z^2) * erfc(-i*z).
im_w_of_x returns Im[w(x)].
The implementation of w_of_z has been entirely rewritten for libcerf-3.6. Depending on the location of z in the complex plane, different expansions are used. For Im z < 0, the reflection formula w(z) = 2 exp(-z^2) - w(-z) is applied, which may amplify the relative error.
If |Im z| is negligible compared to |Re z|, then w(z) is obtained from a first-order expansion around the real axis, where w(x) = exp(-x^2) + i*im_w_of_x(x). Conversely, if |Re z| is negligible compared to |Im z|, then w(z) is obtained from a first-order expansion around the imaginary axis, where w(iy) = erfcx(y). These two cases are treated separately because otherwise the smaller one of the two components of w could be wrong by orders of magnitude, even if the modulus |w| is highly accurate.
For |z| < about 0.23, a Maclaurin series is used.
For about 0.23 < |z| < 7, w(z) is computed from a Taylor expansion of order up to 23 around the center of one of 54 polyominoes that cover this region. The polyominoes, expansion centers and coefficients have been generated by the code generator ttapp (https://jugit.fz-juelich.de/mlz/app/ttapp) for a target relative error of 3 eps = 3*2^-53 in the modulus |w|.
For |z| > 7, an asymptotic expansion is used.
The implementation of im_w_of_x is self-contained. It has been entirely rewritten for libcerf-2.5.
Errors are given in units of eps = 2^-53 = 1.1e-16, as relative deviations of the modulus |w| from high-precision reference values.
For Im z >= 0 and |z| < 7, the relative error of w_of_z is bounded by 3 eps + O(eps^2). This is a proven bound: the Maclaurin series and the tilewise Taylor expansions were generated for it by ttapp, whose error analysis (Wuttke, see REFERENCES) accounts for the truncation and for every rounding, on any platform with IEEE double arithmetic. Measurements stay well inside it: the largest deviation found is 2.3 eps, at the 37622 generated test points of the Taylor domain 0.23 < |z| < 7, up to 20 per tile, as well as at 10^9 random points in that domain, 10^6 per tile; at 10^6 random points in the Maclaurin disc |z| < 0.23 it is 1.2 eps.
For |z| >= 7, the asymptotic expansion has no proven bound yet. The largest deviation found is 4.4 eps, at 3e9 random points with 7 <= |z| < 14, 10^6 per tile; at 26000 random points with |z| up to 1e6 it is 3.2 eps.
For Im z < 0, the reflection formula w(z) = 2 exp(-z^2) - w(-z) is used. Its argument -z^2 = (y-x)(x+y) - 2ixy is rounded to double before the complex exponential is taken, which perturbs the modulus of the result by up to about 1.5 |x^2-y^2| eps and its phase by up to about 0.5 |2xy| eps. The relative error therefore grows in proportion to |z|^2; it stays below about 1.6 |z|^2 eps + 5 eps, and errors close to this bound do occur, for instance 1040 eps at |z| = 28. This limit is intrinsic to the evaluation of exp(-z^2) in double precision. Moreover, near the zeros of w, the first pair of which lies at z = +-1.9915 - 1.3548i, the two terms of the reflection formula cancel, and the relative error is amplified without bound.
The parts Re w and Im w, returned by re_w_of_z and im_w_of_z, are not controlled separately. Where one part is much smaller than the other, its relative error exceeds that of |w| by about the ratio |w|/|Re w| or |w|/|Im w|. This matters above all near the real axis: for y << 1 and x >~ 3, Re w is smaller than |w| by a factor of about y/x, so that the relative error of Re w reaches the order of (x/y) eps; see voigt(3) for the consequences. Within the strips |Im z| < 1e-9 |Re z| and |Re z| < 1e-9 |Im z|, first-order expansions around the axes are used, and the smaller part is accurate to a few |z|^2 eps: up to 2.5 |z|^2 eps were found near the real axis, and up to 4.3 |z|^2 eps near the imaginary axis.
For 0.5 <= |x| < 12, im_w_of_x uses a piecewise Chebyshev approximation generated by ppapp, whose error analysis (Wuttke and Kleinsorge, see REFERENCES) bounds the relative error by 2.2 eps; the largest deviation found at the 432 generated test points in this range is 1.6 eps. For |x| < 0.5 (Maclaurin series) and |x| >= 12 (asymptotic expansion) no theoretical bound exists yet; the relative error was found below 2.7 eps at 816 generated test points for 0.125 <= |x| < 48 and at 30000 random points with |x| up to 1e50.
Milton Abramowitz and Irene M. Stegun, "Handbook of Mathematical Functions", National Bureau of Standards (1964): Formula (7.1.3) introduces the nameless function w(z).
Joachim Wuttke, "Code generation for computing an analytic function with near machine precision on square tiles, with application to the Faddeeva function" (unpublished manuscript, available upon request) explains the tilewise Taylor approximation used for intermediate |z|.
Joachim Wuttke and Alexander Kleinsorge, "Code generation for piecewise Chebyshev approximation", ACM Trans. Math. Softw. 52, 13 (2026), https://doi.org/10.1145/3805698, explains the piecewise polynomial approximation used by im_w_of_x for intermediate values of |x|.
This function is used to compute several other complex error functions: dawson(3), voigt(3), cerf(3), erfcx(3), erfi(3).
Homepage: https://jugit.fz-juelich.de/mlz/lib/cerf
Steven G. Johnson, Massachusetts Institute of Technology, researched the numerics, and implemented the original Faddeeva function.
Joachim Wuttke, Forschungszentrum Juelich, reorganized the code into a library, reimplemented w_of_z using tilewise Taylor approximation and im_w_of_x using piecewise Chebyshev approximation, and wrote this man page.
Please report bugs to the maintainer:
Joachim Wuttke <j.wuttke@fz-juelich.de>
Copyright (c) 2012 Massachusetts Institute of Technology
Copyright (c) 2013-26 Forschungszentrum Juelich GmbH
Software: MIT License.
This documentation: Creative Commons Attribution Share Alike.