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.

The Faddeeva function

$$ 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?


Posted 2017-07-31T12:48:18.843

Reputation: 11 126

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 – 2017-07-31T12:53:35.340

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 – 2017-07-31T13:14:10.717

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.'s ennui – 2017-07-31T13:56:51.833

@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.'s ennui – 2017-07-31T13:58:43.160

@J.M. I'm certainly interested. I've just managed to write up my self-answer before I have to leave work, take a look. – LLlAMnYP – 2017-07-31T14:22:25.563

@LLlAMnYP Sorry for joining the party so late, but these days I want to be sure that someone is really dedicated to a topic before I give an answer that needs more than 10 minutes of my time. I'm providing a fast C-implementation in a moment. – halirutan – 2017-08-18T01:04:27.517

@halirutan thanks, I appreciate any contribution. I'm on vacation till next month, but I'll try to give the topic as much attention as possible. – LLlAMnYP – 2017-08-18T17:55:30.463

I'm looking into my stuff in more detail. Do you need a solution that works to arbitrary precision (slightly hard), or are you fine with a machine precision approximation (easier)? – J. M.'s ennui – 2017-08-20T01:30:44.057

1@J.M arbitrary precision would be extremely interesting from an educational point of view, although for practical purposes it is not necessary. – LLlAMnYP – 2017-08-21T16:32:13.187




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


After that load the package with


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";
FaddeevaWIm::usage = "FaddeevaW[z]";
FaddeevaErfcx::usage = "FaddeevaErfcx[z]";
FaddeevaErfcxRe::usage = "FaddeevaErfcxRe[x]";
FaddeevaErf::usage = "FaddeevaErf[z]";
FaddeevaErfRe::usage = "FaddeevaErfRe[x]";
FaddeevaErfi::usage = "FaddeevaErfi[z]";
FaddeevaErfiRe::usage = "FaddeevaErfi[x]";
FaddeevaErfc::usage = "FaddeevaErfc[z]";
FaddeevaErfcRe::usage = "FaddeevaErfc[x]";
FaddeevaDawson::usage = "FaddeevaDawson[z]";
FaddeevaDawsonRe::usage = "FaddeevaDawsonRe[x]";

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

Let's start with the first one.

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[0]);
    MArgument_setComplex(res, toMComplex(Faddeeva_w(fromMComplex(z), 0)));
    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:

(* 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},
  FaddeevaC = Compile[{{z, _Complex, 0}},
    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:


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]]
Total@Abs[w[small] - FaddeevaC[small]]
Total@Abs[FaddeevaC[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}]]

 measureSerially /@ {w, Faddeeva, FaddeevaC},
 ChartLabels -> {"Mathematica", "LLIAM", "Halirutan"},
 PlotLabel -> "Runtime in seconds for 10^4 values"

Mathematica graphics

It seems we can take the Mathematica implementation out of the equation and only compare yours and mine to see the improvement.

Mathematica graphics

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]]

times = measureParallel /@ {Faddeeva, FaddeevaC}
 ChartLabels -> {"LLIAM", "Halirutan"},
 PlotLabel -> "Parallel Runtime in seconds for 10^7 values"

Mathematica graphics

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

{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.


Posted 2017-07-31T12:48:18.843

Reputation: 109 574

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 – 2017-08-19T10:06:55.073

@LLlAMnYP I appended a small accuracy comparison at the end of my post. – halirutan – 2017-08-20T03:42:20.517

This w[z_] := Exp[-z^2] is surely a typo in the accuracy test? – LLlAMnYP – 2017-08-21T16:30:53.977

@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 – 2017-08-21T19:18:26.170

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 – 2017-09-06T11:14:30.567

The following test works well to find the precision: FadExact[z_] := Exp[-z^2] Erfc[-I z]; exact = Rationalize[RandomComplex[{-5 - 5 I, 5 + 5 I}, 10^4], 0]; {outHal, outLLL, outExact} = {FaddeevaW[exact], Faddeeva[exact], N[FadExact[exact], 30]}; outHal/outExact//Abs//Log//Max and outLLL/outExact//Abs//Log//Max both return something just under 3.*^-14, so the precision is 13-14 digits even for negative imaginary parts. I'm quite convinced by your implementation, so I'm accepting your answer. Thanks a lot! – LLlAMnYP – 2017-09-07T10:51:39.463

@LLlAMnYP I didn't anticipate that such a compilation error will occur. In fact both on Linux and Mac OS, CreateLibrary uses the -O2 switch per default and the -std=c99 option is not necessary. Thanks for accepting the answer. I have written it in the hope that you will consider writing a C-wrapper yourself next time when you a good C-library at hand that does the job. Setting up a first usable library-function did not take more than half an hour. – halirutan – 2017-09-07T13:28:10.103

Thanks for your contribution. I presume, this is to do with the value of $CCompiler which, when everything is set up properly, likely has the necessary switch, but on my system I only have a portable compiler and had to set up $CCompiler manually. I'll certainly use your input in the future as a reference for writing C-wrappers -- the more advanced subjects in MMA are often not so well documented. As an aside, I find it curious, that ~2k lines of C code are so much faster than a relatively small function that is compiled to C anyway. – LLlAMnYP – 2017-09-08T11:14:40.080

I have gotten round to properly learning LibraryLink and of course, listability came to mind. I didn't pay heed to your mention of parallelization until today I came across your answer here. Now above you have faddeevaCCall = LibraryFunctionLoad[... but in the next code blocks you write With[{fc = faddeevaC},... I assume this was just a typo and faddeevaCCall should simply be faddeevaC? Anyway, I'll leave a useful link in the body of the Q

– LLlAMnYP – 2018-05-18T13:27:19.187

@LLlAMnYP It's a typo. It should be With[{fc = faddeevaCCall}, ..]. Please look in the linked package code. It's a tiny 100 lines package and there, it is correct. What I do is wrap the library function inside another Compile to make it parallelized and listable. That's the entire trickery :) – halirutan – 2018-05-18T13:58:15.793

I thought so too. Very neat that the wrapper doesn't cause a slowdown. – LLlAMnYP – 2018-05-18T14:05:31.233



The current state of my code is available at

Minor updates (03.08.17):

  • 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},
          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}] // Developer`ToPackedArray,
        lookupAN = Delete[Table[n/4., {n,-106, 106}],{107}], lookupA2N2},
        lookupEmA2N2[[3]] *= I lookupAN;
        lookupA2N2 = lookupAN^2 + 0.I;
    With[{lookup = lookupEmA2N2//Developer`ToPackedArray, lAN = lookupAN//Developer`ToPackedArray, lA2N2 = lookupA2N2//Developer`ToPackedArray},
        {{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[[2]] *= e2ANx; sums[[3]] *= 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} $$


$$\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.


Posted 2017-07-31T12:48:18.843

Reputation: 11 126

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.'s ennui – 2017-07-31T14:45:52.140

@J.M. They are merged :-) – LLlAMnYP – 2017-07-31T14:58:38.873


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]:

real and imaginary parts of difference

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[10], {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[9]/18]), 
                 ContourShading -> None, ContourStyle -> ColorData[97, 2], 
                 PlotPoints -> 55], AspectRatio -> Automatic]

altitude chart

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)$:

domain coloring of w(z)

J. M.'s ennui

Posted 2017-07-31T12:48:18.843

Reputation: 115 520

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 – 2017-08-24T13:50:19.000

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.'s ennui – 2017-08-24T14:22:11.500


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]:

relative error

The altitude chart, clothoid curve, and domain colored image in my other answer are reproduced by this routine.

J. M.'s ennui

Posted 2017-07-31T12:48:18.843

Reputation: 115 520

I could not make it work in mathematica 11.3, something with With[] syntax :( – Zlelik – 2018-09-30T10:06:05.173

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.'s ennui – 2018-09-30T14:51:21.203

Yes, I can wait. – Zlelik – 2018-09-30T21:11:44.187

@Zlelik, please try the routine now – J. M.'s ennui – 2018-10-01T18:13:37.063

Thanks a lot. Now it works. It 75% faster than standard function. – Zlelik – 2018-10-02T19:54:12.347