# Implementing a compilable Faddeeva function of complex argument

For those, who are looking at Halirutan's answer and thinking "gee, I wish I was that good at LibraryLink, then I could really speed up my code!" I leave here the following links:

I think I mostly got the hang of using LibraryLink in the space of a day with these links, especially the first one.

$$w(z) = \exp(-z^2)(1 - \operatorname{erf}(-iz)) = \exp(-z^2)\operatorname{erfc}(-iz)=\operatorname{erfcx}(-iz)$$

has important applications in the field of optical spectroscopy. It is, in particular, related to the real and imaginary Voigt functions as

$$w(x + i y) = V(x,y) + iL(x,y).$$

The various error functions of complex argument, expressible in Mathematica as DawsonF, Erf, Erfi, and Erfc are not compilable (Erfi, DawsonF) or compilable only for real arguments. My goal is to implement the Faddeeva function as a compiled function of complex argument for the purposes of, e.g., better performance inside Dynamic or for fitting of complex-valued datasets.

As a starting point, I have considered the Faddeeva package which provides an implementation in C/C++ code with wrappers for various languages, among which, unfortunately, Mathematica is not present.

The algorithm used in the linked package is based on the following pre-print:

which works well for relatively small-magnitude values of $$z$$, and a continued fraction approximation for large values. I, unfortunately, do not have access to this article, but the approximation is given by the following Mathematica code:

With[{max = 2*n},
Fold[(#2/2)/(z - #) &, (max + 1)/2/z,
Range[max, 1, -1] ~Join~ {2/Sqrt[Pi]}]] // Together


which for integer n gives:

$$\frac{2 z}{\sqrt{\pi}} \cdot \frac{P(z^2)}{Q(z^2)}$$

where $$P$$ and $$Q$$ are polynomials of degree n and n+1, respectively.

This approximation shows very good convergence already for Abs[z] > 10. Obviously, this is the easy part. The much harder part is constructing a decent algorithm for smaller z.

How can I achieve this? / How can I improve my algorithm?

• For small absolutely z, you might be able to use a Padé expansion to obtain a ratio of polynomials. Have a look at the documentation page of PadeApproximant. – Henrik Schumacher Jul 31 '17 at 12:53
• Another option is to write a wrapper for the c++ code with the help of LTemplate. This is not as complicated as it seems. – Henrik Schumacher Jul 31 '17 at 13:14
• I have a routine, but it needs to be cleaned up a bit. I hope you don't mind a wait of a few days. These two are related questions. – J. M. will be back soon Jul 31 '17 at 13:56
• @Henrik, in the particular case of the error/Voigt/Faddeeva function, one must first factor out the exponential behavior before constructing the Padé approximant; otherwise the approximant is numerically unstable. – J. M. will be back soon Jul 31 '17 at 13:58
• @J.M arbitrary precision would be extremely interesting from an educational point of view, although for practical purposes it is not necessary. – LLlAMnYP Aug 21 '17 at 16:32

# Preface

I've put a FaddeevaM package on my GitHub account that will do all that follows automatically. It is a full package that exposes all (and not only the Faddeeva $w$) functions for Mathematica users.

I have tested the package under Linux (10.3 and 11.1) and Mac OS X (11.1) and it seems to work. I will put further documentation in the GitHub repository. To install it you can download the paclet-file from the releases section and use

PacletInstall["path/to/downloaded/file.paclet"]


After that load the package with

<<FaddeevaM


Here comes the crucial part. I provided a compiled library for Linux (and I can include the one for OS X later). If the library is not available for your system, my package tries to compile it from source. This only works if you have a usable C-compiler installed and you can successfully evaluate the example from the CreateLibrary documentation.

The package will tell you when it tries to rebuild the code and show you the compilation command. If no further error-messages appear, everything is set up. This needs only to be done once.

The package provides the following functions:

FaddeevaW::usage = "FaddeevaW[z] calculated Exp[-z^2] (1 - Erf[-I*z]) for complex z";


Create an issue on GitHub if something goes wrong.

Here we go: Grab the C-files and header from the Faddeeva Package and put it in a directory on your disk. What you need is an additional c-file that acts as wrapper so that LibraryLink can call the functions from the package.

Create your wrapper MMAFaddeeva.c or you use mine. The basic ingredients are 2 things:

• convert from Mathematica C complex type to the standard complex type and back
• call the package function and return the result

inline complex double fromMComplex(mcomplex z) {
return mcreal(z)+I*mcimag(z);
}

inline mcomplex toMComplex(complex double z) {
mcomplex result = {creal(z), cimag(z)};
return result;
}


The wrapper function that is then loaded in Mathematica is easy as well

DLLEXPORT int m_Faddeeva_w(WolframLibraryData libData, mint argc, MArgument *args, MArgument res) {
mcomplex z = MArgument_getComplex(args);
return LIBRARY_NO_ERROR;
}


Now, you can create a library from the code directly in Mathematica, where files are the paths to your c-files (check with FileExistQ).

lib = CreateLibrary[files, "libFaddeeva"]


If the library was created successfully, you can load the function

faddeevaCCall = LibraryFunctionLoad[lib, "m_Faddeeva_w", {_Complex}, _Complex];


This is directly usable:

faddeevaCCall
(* 0.367879 + 0.607158 I *)


but we probably want this function to be Listable so that it can be called in parallel. To do this, we wrap it with a Compile that does the work for us:

With[{fc = faddeevaCCall},
fc[z],
RuntimeAttributes -> {Listable},
Parallelization -> True
]
];


Now we can run some tests. First, I'm defining a high-level Mathematica function of w and I grab your implementation from GitHub:

Get["https://raw.githubusercontent.com/lllamnyp/Faddeeva/master/Faddeeva.m"];

SetAttributes[w, {Listable}];
w[z_] := Exp[-z^2] (1 - Erf[-I*z]);

small = RandomComplex[.3 {-15 - 15 I, 15 + 15 I}, 10^4];


Let's compare the total error for the 10^4 random numbers.

Total@Abs[w[small] - Faddeeva[small]]
(* 0.0000357799 *)
(* 0.000036873 *)
(* 4.59592*10^-6 *)


It seems your implementation gives similar results as the C-library. Let's look at the speed when we call the function serially:

measureSerially[f_] := First@AbsoluteTiming[Do[f[z], {z, small}]]

BarChart[
ChartLabels -> {"Mathematica", "LLIAM", "Halirutan"},
PlotLabel -> "Runtime in seconds for 10^4 values"
] It seems we can take the Mathematica implementation out of the equation and only compare yours and mine to see the improvement. On my machine, this is a speedup of about 30x. Let's do some real work and use a large sample calling both our functions in parallel

large = RandomComplex[.3 {-15 - 15 I, 15 + 15 I}, 10^7];

measureParallel[f_] := First@AbsoluteTiming[f[large]]

BarChart[
times,
ChartLabels -> {"LLIAM", "Halirutan"},
PlotLabel -> "Parallel Runtime in seconds for 10^7 values"
] Looks equally, but on my machine, the speedup is now over 42x.

# Accuracy comparison

For finding out the accuracy in comparison to LLlAMnYP's approach, let us define some rational values and calculate them with high precision. For this, we use the standard definition of Faddeeva $w$:

data = Table[x, {x, -10, 10, 1/100}];
SetAttributes[w, {Listable}];
w[z_] := Exp[-z^2] (1 - Erf[-I*z]);


First, median and maximum absolute error of LLlAMnYP's approach

Get["https://raw.githubusercontent.com/lllamnyp/Faddeeva/master/Faddeeva.m"];
{Median[#], Max[#]} &@Abs[N[w[data], 40] - Faddeeva[data]]
(* {2.498*10^-16, 1.56486*10^-13} *)


Now the functions from the Faddeeva package

{Median[#], Max[#]} &@Abs[N[w[data], 40] - FaddeevaW[data]]
(* {1.38778*10^-17, 2.23773*10^-16} *)


It seems we gained about one order of magnitude in precision in the median and several more in the maximum error.

• While I'm not able to check in depth right now, there's a minor comment about testing for precision and accuracy. Try evaluating the top level w with rational arguments and arbitrary precision,such as N[..., 30]. Otherwise mma's machine precision calculations are way off. – LLlAMnYP Aug 19 '17 at 10:06
• @LLlAMnYP I appended a small accuracy comparison at the end of my post. – halirutan Aug 20 '17 at 3:42
• This w[z_] := Exp[-z^2] is surely a typo in the accuracy test? – LLlAMnYP Aug 21 '17 at 16:30
• @LLlAMnYP Yes, this somehow slipped in. I re-ran the tests with the correct definition and got the same results as shown above. It must have been a copy and paste mistake I made. – halirutan Aug 21 '17 at 19:18
• I'm back from vacation and have gotten around to testing this properly. The library wouldn't build until I added "CompileOptions" -> "-O2 -std=c99" to $CCompiler$, but this was through some trial and error. Could you perhaps touch a bit on the subject? I'm not sure if I'm doing things in the cleanest fashion. – LLlAMnYP Sep 6 '17 at 11:14

TL;DR

The current state of my code is available at

• I replaced the Total[Divide[..., ...]] in the definition of Faddeeva for Abs[z] < 10 with a string of two Dot products, it is now slightly faster.

• When passing CompilationTarget -> "C" to Faddeeva I kept getting a bunch of ::cmperr warnings, as a result a LibraryFunction was not being generated. After removing RuntimeOptions -> "Speed" it now works properly and can be inlined in further compiled functions.

• I disabled symbolic evaluation. Now things like Plot[Faddeeva[x] // ReIm, {x, -5, 5}, Evaluated -> True work properly.

As do the authors of the aforementioned Faddeeva package, I begin by implementing the scaled real and imaginary error functions of real arguments:

$$\mathrm{ErfcxRe}[x] = \mathrm{Exp}[x^2]\mathrm{Erfc}[x],$$

$$\mathrm{ErfcxIm}[x] = \mathrm{Exp}[-x^2]\mathrm{Erfi}[x].$$

Later on I realized that the former is compilable and seems to be faster anyhow, so I ended up not using it in the implmentation of Faddeeva, but it is nonetheless present in the package. It is implemented in a manner quite similar to the implementation of ErfcxIm, as I detail below.

ErfcxIm = With[{lookupErfcxIm = lookupErfcxIm},
Compile[{{x, _Real}},
Block[{res = 0., xAbs = Abs[x]},
If[xAbs >= 5.*^7, Return[Divide[0.5641895835477563,x]]];
If[xAbs > 48.9, With[{x2=x*x},
Return[Divide[
0.5641895835477563 x (-558+x2 (740+x2 (-216+16 x2))),
105+x2 (-840+x2 (840+x2 (-224+16 x2)))]]]];
If[x == 0., Return[res]];
res = With[{lookupTable = lookupErfcxIm, y = Divide[1, 1 + xAbs]},
With[{n = Subtract[Floor[100 y], 1]},
With[{yOff = y - 0.015 - 0.01 n},
Fold[# yOff + #2 &, Reverse[lookupTable[[n]] ] ]
]]];
res * Sign[x]
],
RuntimeAttributes->{Listable}, RuntimeOptions->"Speed",
Parallelization->True, CompilationTarget->"C"]
];


0.5641895835477563 is 1/Sqrt[Pi]. For very large arguments (>5e7) a 1/x dependence is within machine precision, for smaller ones the continued fraction approximation works well, for the yet smaller arguments the function is reparametrized as follows (same approach as in the linked Faddeeva C-package): y -> 1/(1+x), which is subdivided into several (50) intervals and each is fitted with a 7th degree polynomial of y such that the adjacent intervals have matching derivatives up to 3rd order. The polynomials' are taken in HornerForm and their coefficients are stored in the lookupTable, a {50, 8} matrix. The polynomial is reconstructed by the Fold statement. Using MMA's arbitrary precision arithmetic I have verified that this gives results very close to machine precision.

This function is useful for one special case. It is used for calculating the Faddeeva function of purely real arguments (in optical spectroscopy this corresponds to a Gaussian absorption line).

The general form of the Faddeeva function for smallish complex arguments is, as given in the linked pre-print, the following (I have factored out the $\exp(-x^2)$ from the sums in the following):

$$Re[w(z)] = e^{-x^2} \left\{\mathrm{ErfcxRe}(y) \cos(2xy) + \\ 2 a \left[ x \sin(xy) \mathrm{sinc}(xy) -y \left(\cos(2xy)\Sigma_1 - \Sigma_2/2 - \Sigma_3/2\right)\right]/\pi\right\}$$

$$Im[w(z)] = e^{-x^2} \left\{-\mathrm{ErfcxRe}(y) \sin(2xy) + \\ 2 a \left[ x \mathrm{sinc}(2xy) + y \sin(2xy)\Sigma_1 - \Sigma_4/2 + \Sigma_5/2\right]/\pi\right\}$$

the $\Sigma_i$ are

$$\left\{\Sigma_1, \Sigma_2, \Sigma_3, \Sigma_4, \Sigma_5\right\} = \\ \sum_{n=1}^\infty \frac{\exp(-a^2n^2)}{a^2n^2+y^2}\left\{1, e^{-2anx}, e^{2anx}, an e^{-2anx}, an e^{2anx} \right\}$$

I slightly simplify the sums by redefining:

$$\Sigma_1 = \frac{1}{2}\sum_{n=-\infty}^\infty \frac{\exp(-a^2n^2)}{a^2n^2+y^2} = \Sigma$$

$$\Sigma_2 + \Sigma_3 = \sum_{n=-\infty}^\infty \frac{\exp(2anx-a^2n^2)}{a^2n^2+y^2} = \Sigma_{23}$$

$$\Sigma_5 - \Sigma_4 = \sum_{n=-\infty}^\infty \frac{an \exp(2anx-a^2n^2)}{a^2n^2+y^2} = \Sigma_{45}$$

In all the above sums n==0 is excluded.

There are several corresponding terms in the real and imaginary parts which can be nicely converted to exponential form:

$$w(z) = e^{-x^2} \left\{\mathrm{ErfcxRe}(y)e^{-2ixy} + a\left[2 ix\, \mathrm{sinc}(xy)e^{-ixy} - y e^{-2ixy} \Sigma + y \Sigma_{23} + i\Sigma_{45}\right]/\pi\right\}$$

$a$ is a parameter used in an approximation for $e^{t^2}$ to provide an analytical solution for the integral involving $e^{t^2}$. The pre-print shows nicely that a<.5 is by far sufficient to achieve machine precision. In my implementation I use a==.25.

In order to speed up computation, I precompute tables of $\exp(-a^2n^2)$, $an$, $a^2n^2$ which are injected into the body of the compiled function with a With. I take -106<=n<=106. This almost reaches $MinMachineNumber for Exp[-a^2 n^2]. With all these definitions out of the way, the definition of the Faddeeva function is as follows: Faddeeva = Block[{lookupEmA2N2 = Table[Reverse@#~Join~#&[Table[Exp[-n^2/16.] + 0. I,{n,106}]],{3}] // DeveloperToPackedArray, lookupAN = Delete[Table[n/4., {n,-106, 106}],{107}], lookupA2N2}, lookupEmA2N2[] *= I lookupAN; lookupA2N2 = lookupAN^2 + 0.I; With[{lookup = lookupEmA2N2//DeveloperToPackedArray, lAN = lookupAN//DeveloperToPackedArray, lA2N2 = lookupA2N2//DeveloperToPackedArray}, Compile[ {{z, _Complex}},With[{x = Re[z], y = Im[z]},With[{ere = Exp[y*y]Erfc[y], mxx = Minus[x*x], mxy = Minus[x*y]}, If[z == 0. + 0. I, Return[1. + 0. I]]; If[x == 0., Return[0. I + ere]]; If[y == 0., Return[Exp[mxx] + I ErfcxIm[x]]]; If[Abs[z] < 10., Block[{sums = lookup, e2ANx = Exp[2 lAN x]}, sums[] *= e2ANx; sums[] *= e2ANx; Return[Exp[mxx] * (ere Exp[2I mxy] + 0.07957747154594767 * (2 I x Sinc[mxy] Exp[I mxy] + {Minus[y Exp[2 I mxy]], y, 1.}.sums.Divide[1,(lA2N2 + y*y)])) ] ] ]; With[{zz=z*z},0.5641895835477563 I Divide[z (-558+zz (740+zz (-216+16 zz))),105+zz (-840+zz (840+zz (-224+16 zz)))]]]], RuntimeAttributes->{Listable}, RuntimeOptions->{"EvaluateSymbolically"->False}, Parallelization->True, CompilationOptions->{"InlineExternalDefinitions" -> True, "InlineCompiledFunctions"->False}, CompilationTarget-> "C" ] ] ]  Within the body of the last Return statement for Abs[z]<10 sums evaluates to a table of the form $$sums = \left\{ \exp(-a^2n^2), \exp(2anx-a^2n^2), ian \exp(2anx-a^2n^2) \right\}_{n=-106..106}$$ so $$\left\{-ye^{-2ixy},y,1\right\}.sums.\left( 1/\left\{a^2n^2+y^2\right\}\right)$$ gives the last 3 addends of$w(z)$in the square braces. A quick performance check: test = RandomComplex[.3 {-15 - 15 I, 15 + 15 I}, 10^4]; Faddeeva[test] // AbsoluteTiming // First (* 0.149745 *) Exp[-test^2] Erfc[-I test] // AbsoluteTiming // First (* 2.88898 *)  20x speed increase. Could this be better? After all, .1ms per complex number in the difficult range is not that great. A precision test: N[Table[Exp[-(x + I y)^2] Erfc[-I (x + I y)], {x, 0, 5, 1/20}, {y, 0, 5, 1/20}], 30]/ Table[Faddeeva[x + I y], {x, 0, 5, 1/20}, {y, 0, 5, 1/20}] // Log // Abs // Max (* 6.63025*10^-13 *)  Roughly 12-13 digits of precision at worst. Certainly good enough for experimental data. • Don't have time to read it thoroughly, but: since Mathematica supports complex inputs and outputs, you might consider merging the routines for the real and imaginary parts so that genuine complex output is returned, especially since the computations for both seem to have common quantities. – J. M. will be back soon Jul 31 '17 at 14:45 • @J.M. They are merged :-) – LLlAMnYP Jul 31 '17 at 14:58 Here is a compiled implementation of the Faddeeva function$w(z)$, based on the algorithm of Abrarov and Quine (see this and this for reference): faddeeva = With[{sp = N[Sqrt[π]], τ = 18, nn = Ceiling[18^2/(2 π)]}, With[{cof = Table[2 Exp[-(π k/τ)^2]/(sp τ (sp k)^2), {k, nn}], d = N[π] Range[nn], s = PadRight[{}, nn, {-1., 1.}]}, Compile[{{z, _Complex}}, Module[{zz = z, f = 1. + 0. I, l = Im[z] < 0., e, w}, If[z != 0., If[l, zz = -z]; w = τ zz; e = Exp[I w]; f = I ((1 - e)/w + (τ w/sp) cof.((s e - 1)/(1 - (w/d)^2))); If[l, f = 2 Exp[-z^2] - f]]; f], RuntimeAttributes -> {Listable}]]];  The routine works in principle for all complex arguments; in practice, the computation is less numerically stable for$\Im z<0$. Nevertheless, a lot of the applications of$w(z)$have$\Im z\ge 0$(e.g. the Voigt function), so this may not be much of a disadvantage. faddeeva is nearly a hundred times as efficient as an implementation in terms of Erfc[], without even compiling to C: RepeatedTiming[faddeeva[Table[x + I y, {x, -8, 8, 1/20}, {y, -8, 8, 1/20}]];] {0.27, Null} fw[z_] := Exp[-z^2] Erfc[-I z]; RepeatedTiming[fw[N @ Table[x + I y, {x, -8, 8, 1/20}, {y, -8, 8, 1/20}]];] {26.6, Null}  Here is a plot showing the real and imaginary parts of faddeeva[z] - Exp[-z^2] Erfc[-I z]: Here is a colorized version of the "altitude chart" of$w(z)$as seen in Abaramowitz and Stegun: Show[ContourPlot[Abs[faddeeva[x + I y]], {x, 0, 3 + 1/5}, {y, -3, 3}, Contours -> Join[Subdivide, {2, 3, 4, 5, 10, 100}], ContourShading -> None, ContourStyle -> ColorData[97, 1], PlotPoints -> 55], ContourPlot[Arg[faddeeva[x + I y]], {x, 0, 3 + 1/5}, {y, -3, 3}, Contours -> (π Union[Subdivide[-1, 1, 6], Range/18]), ContourShading -> None, ContourStyle -> ColorData[97, 2], PlotPoints -> 55], AspectRatio -> Automatic] Here is a plot of a clothoid: ParametricPlot[ReIm[Exp[I t^2] faddeeva[t Exp[I π/4]]], {t, -7, 7}, Axes -> None, Frame -> True] Finally, here is a domain-colored image of$w(z)$: • Thanks for this contribution, I'll clearly have to spend some time studying the algorithm. IMO, both your and halirutan's solutions are equally worthy of the bounty, but after awarding the first one, I'm somehow no longer able to offer you the 250 rep. – LLlAMnYP Aug 24 '17 at 13:50 • It's fine; I think his answer deserves a bounty. I've added more examples in the meantime. (Also, I still have to work on the arbitrary precision version, tho that one is not really compilable anymore.) – J. M. will be back soon Aug 24 '17 at 14:22 Prompted by this related question, I decided to revisit the numerical evaluation of the Faddeeva function. As it turns out, Abrarov and Quine, the persons who came up with the algorithm in my other answer, have recently developed a new method with better coverage. With a few modifications, here is an adaptation of their code to Mathematica: faddeeva = With[{n = 23, σ = 1.5, h = N[3/23/π], f = N[2/Sqrt[π]], eps =$MachineEpsilon},
With[{alf = 8 π h Table[k Exp[-(2 π h k)^2] Sin[4 π h k σ], {k, n}],
bet = 2 Table[Exp[-(2 π h k)^2] Cos[4 π h k σ], {k, n}],
cn = (2 π h Range[n])^2, he = h Exp[σ^2]},
Compile[{{z, _Complex}},
Module[{zf = z, c, d, g, k, u, v, w, zs},
If[Im[zf] < 0, zf = -zf];
If[Abs[zf] <= 7,
If[Im[zf] < 0.1,
zs = zf^2 + σ^2;
w = Exp[-zf^2] + 2 I he zf
(1/zs + (σ alf + bet (zs - cn)).
(1/(4 σ^2 cn + (cn - zs)^2))),
zs = zf + I σ;
w = 2 I he/zs + he (alf - 2 I zs bet).
(1/(cn - zs^2))],
w = v = 2. zf^2 - 1.;
c = w; d = 0. I; k = 1;
While[u = k (k + 1); v -= 4.;
c = v - u/c; d = 1/(v - u d);
w *= (g = c d); k += 2;
Abs[g - 1] > eps];
w = I zf f/w];
If[Im[z] < 0, w = 2 Exp[-z^2] - w]; w],
RuntimeAttributes -> {Listable}]]];


This implementation is a little different from the one by Abrarov and Quine. As with their method, for $$|z|\le7$$, the two different approximations they derived are used; outside that disc, I use the Lentz-Thompson-Barnett algorithm (instead of backward evaluation in the original) on the even part of the Laplace continued fraction used in the original.

Tests with RepeatedTiming[], as was done in my other answer, show a significant speedup.

Here is a plot of the real and imaginary parts of the relative error against fw[z_] := Exp[-z^2] Erfc[-I z]: The altitude chart, clothoid curve, and domain colored image in my other answer are reproduced by this routine.

• I could not make it work in mathematica 11.3, something with With[] syntax :( – Zlelik Sep 30 '18 at 10:06
• Ah, the With[] there is undocumented syntax, sorry. :( If you can wait, I will edit this later to use the more conventional syntax. – J. M. will be back soon Sep 30 '18 at 14:51
• Yes, I can wait. – Zlelik Sep 30 '18 at 21:11
• @Zlelik, please try the routine now – J. M. will be back soon Oct 1 '18 at 18:13
• Thanks a lot. Now it works. It 75% faster than standard function. – Zlelik Oct 2 '18 at 19:54