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 );
The data type double complex is defined in the header <complex.h>, which C99 introduced. Since C11 it is optional: an implementation may define __STDC_NO_COMPLEX__ and then provide neither the header nor the type, and Microsoft's C compiler does not support the arithmetic operators for it. As a fallback, use the C++ variant of this library, libcerfcpp, in which double complex is replaced by std::complex<double> from the header <complex>.
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 with up to 23 terms is used. It is summed by a real three-term recurrence (the device of Goertzel and Clenshaw), which needs no complex division and yields Re w and Im w free of cancellation.
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 rounding error of the asymptotic expansion is bounded by 5.1 eps for the modulus and 5.4 eps for the components, again a proven bound (Wuttke, see REFERENCES). At 200000 random points with |z| up to 1e6 the largest deviation of the modulus found is 3.7 eps, and Re w and Im w taken separately deviate by up to 4.5 eps. Note that these two figures bound the rounding error alone. Adding the truncation error, which is below 0.1 eps of |w|, gives a bound for the modulus, but not for the components: for that, a lower bound on each component separately would be needed, and only the one for Re w is proven so far.
For Im z < 0, the reflection formula w(z) = 2 exp(-z^2) - w(-z) is used. Its exponent -z^2 = (y-x)(x+y) - 2ixy is carried as an unevaluated sum of two doubles, so that the factor is obtained to about 4 eps whatever |z| (Wuttke, see REFERENCES). At 20000 random points where the exponential neither over- nor underflows, |z| reaching 1e8, the largest deviation of the modulus found is 5.4 eps. Before libcerf-3.7 the exponent was rounded into a single double, and the relative error grew in proportion to |z|^2, reaching 1040 eps at |z| = 28 and leaving no correct digit at all where x^2-y^2 stays small up to |z| = 1e8. 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|, which is unbounded. This matters above all near the real axis, where Re w is smaller than |w| by a factor of about y/x; see voigt(3) for the consequences.
Three regions have to be distinguished, and only the first two are controlled.
For |z| < 7 and |y| < 1e-2 |x|, Re w is computed from an expansion about the real axis rather than from the tilewise Taylor approximation, whose bound covers the modulus alone. Terms of that expansion shrink by (y/x)^2, and five of them are carried, so that the truncation is far below the last bit. Its first coefficient is the derivative of Im w, which has its own approximation and is accurate to under 4 eps; formed instead as 2/sqrt(pi) - 2 x Im w(x) it would cancel to relative size 1/(2 x^2) and cost 2 x^2 eps. The relative error of Re w stays below 16 eps, measured over a scan of 9800 points, against 6e8 eps before this expansion was introduced. Im w, and with it the modulus, still comes from the Taylor tiles and keeps the bound of 3 eps stated above.
Within the strip |Im z| < 1e-9 |Re z| a first-order expansion about the real axis is used, with the same derivative as its coefficient; the relative error of Re w stays within a few eps there. Within the mirror strip |Re z| < 1e-9 |Im z| the corresponding expansion has the derivative of erfcx as its coefficient, which likewise has an approximation of its own, accurate to under 3.8 eps; formed instead as 1/sqrt(pi) - y erfcx(y) it would cancel to relative size 1/(2 y^2) and cost 2 y^2 eps. The relative error of Im w stays below 1.2 eps there, measured at |z| from 2 to 1e5 and |Re z| / |Im z| from 1e-10 to 1e-20. Outside that strip Im w keeps its relative accuracy, because the tables of the tilewise Taylor approximation place their expansion centres on the imaginary axis, which preserves the symmetry w(z*) = w(-z)*.
For |z| >= 7 and |y| << |x| the asymptotic expansion is used. Being a power series in 1/z it cannot supply the term exp(-x^2), which is all of Re w on the axis, so the relative error of Re w grows once y falls below about sqrt(pi) x^2 exp(-x^2) / eps. That threshold exceeds the strip limit 1e-9 x only for 7 <= x <~ 7.75; there, up to 4e4 eps were found, at x = 7, y = 1e-8. Elsewhere the strip already covers the affected region.
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, "libcerf, complex error function and related functions reimplemented with relative accuracy guarantees" (unpublished manuscript, available upon request) derives the error bound of the asymptotic expansion used for large |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.