evaluate Faddeeva function for complex argument
- Author
- Manuel Schiller manue.nosp@m.l.sc.nosp@m.hille.nosp@m.r@ni.nosp@m.khef..nosp@m.nl
- Date
- 2013-02-21
Calculate the value of the Faddeeva function w(z) = \exp(-z^2) \mathrm{erfc}(-i z).
The method described in
S.M. Abrarov, B.M. Quine: "Efficient algotithmic implementation of
Voigt/complex error function based on exponential series approximation" published in Applied Mathematics and Computation 218 (2011) 1894-1902 doi:10.1016/j.amc.2011.06.072
is used. At the heart of the method (equation (14) of the paper) is the following Fourier series based approximation:
w(z) \approx \frac{i}{2\sqrt{\pi}}\left( \sum^N_{n=0} a_n \tau_m\left( \frac{1-e^{i(n\pi+\tau_m z)}}{n\pi + \tau_m z} - \frac{1-e^{i(-n\pi+\tau_m z)}}{n\pi - \tau_m z} \right) - a_0 \frac{1-e^{i \tau_m z}}{z} \right)
The coefficients a_b are given by:
a_n=\frac{2\sqrt{\pi}}{\tau_m} \exp\left(-\frac{n^2\pi^2}{\tau_m^2}\right)
To achieve machine accuracy in double precision floating point arithmetic for most of the upper half of the complex plane, chose t_m=12 and N=23 as is done in the paper.
There are two complications: For Im(z) negative, the exponent in the equation above becomes so large that the roundoff in the rest of the calculation is amplified enough that the result cannot be trusted. Therefore, for Im(z) < 0, the symmetry of the erfc function under the transformation z --> -z is used to avoid accuracy issues for Im(z) < 0 by formulating the problem such that the calculation can be done for Im(z) > 0 where the accuracy of the method is fine, and some postprocessing then yields the desired final result.
Second, the denominators in the equation above become singular at z = n * pi / 12 (for 0 <= n < 24). In a tiny disc around these points, Taylor expansions are used to overcome that difficulty.
This routine precomputes everything it can, and tries to write out complex operations to minimise subroutine calls, e.g. for the multiplication of complex numbers.
In the square -8 <= Re(z) <= 8, -8 <= Im(z) <= 8, the routine is accurate to better than 4e-13 relative, the average relative error is better than 7e-16. On a modern x86_64 machine, the routine is roughly three times as fast than the old CERNLIB implementation and offers better accuracy.
For large |z|, the familiar continued fraction approximation
w(z)=\frac{-iz/\sqrt{\pi}}{-z^2+\frac{1/2}{1+\frac{2/2}{-z^2 + \frac{3/2}{1+\frac{4/2}{-z^2+\frac{5/2}{1+\frac{6/2}{-z^2+\frac{7/2 }{1+\frac{8/2}{-z^2+\frac{9/2}{1+\ldots}}}}}}}}}}
is used, truncated at the ellipsis ("...") in the formula; for |z| > 12, Im(z)>0 it will give full double precision at a smaller computational cost than the method described above. (For |z|>12, Im(z)<0, the symmetry property w(x-iy)=2e^{-(x+iy)^2-w(x+iy)} is used.
Definition at line 542 of file RooMath.cxx.