# erf(a+ib) error function separate into real and imaginary part

Is there an easy way to separate erf(a+ib) into real and imaginary part?

I'm not sure if you are interested in an analytical answer or a computational answer; these are two different things. The analytical answer is...not really, unless you consider GEdgar's answer useful. (And one might.)

The computational answer is a resounding yes. A result found in Abramowitz & Stegun claims the following:

$$\operatorname*{erf}(x+i y) = \operatorname*{erf}{x} + \frac{e^{-x^2}}{2 \pi x} [(1-\cos{2 x y})+i \sin{2 x y}]\\ + \frac{2}{\pi} e^{-x^2} \sum_{k=1}^{\infty} \frac{e^{-k^2/4}}{k^2+4 x^2}[f_k(x,y)+i g_k(x,y)] + \epsilon(x,y)$$

where

$$f_k(x,y) = 2 x (1-\cos{2 x y} \cosh{ k y}) + k\sin{2 x y} \sinh{k y}$$ $$g_k(x,y) = 2 x \sin{2 x y} \cosh{k y} + k\cos{2 x y} \sinh{k y}$$

Then

$$\left |\epsilon(x,y) \right | \le 10^{-16} |\operatorname*{erf}{(x+i y)}|$$

This accuracy is valid for all $x$ and $y$, i.e., the complex plane.

I will present a derivation of this result to show you where the error term comes from. Consider the definition of the error function in the complex plane:

$$\operatorname*{erf}{z} = \frac{2}{\sqrt{\pi}} \int_{\Gamma} d\zeta \, e^{-\zeta^2}$$

where $\Gamma$ is any path in the complex plane from $\zeta = 0$ to $\zeta=z$. Consider, then, the special case where $\Gamma$ is the path that runs from $0$ to $x$ along the real axis, then from $x$ to $z=x+i y$ parallel to the imaginary axis. Seen this way, the error function of a complex number is equal to

$$\operatorname*{erf}{(x+i y)} = \operatorname*{erf}{x} + i \frac{2}{\sqrt{\pi}} e^{-x^2} \int_0^y du \, e^{u^2} \cos{2 x u} \\ + \frac{2}{\sqrt{\pi}} e^{-x^2} \int_0^y du \, e^{u^2} \sin{2 x u}$$

In what follows, we will find an appropriate approximation to these integrals. To do this, we take a detour through some Fourier theory.

Consider a function $\phi(t)$ that has a Fourier transform

$$\Phi(\xi) = \int_{-\infty}^{\infty} dt \, \phi(t) \, e^{-i 2 \pi \xi t}$$

We begin with a form of the Poisson sum formula:

$$\sum_{n=-\infty}^{\infty} \phi(n+t) = \sum_{n=-\infty}^{\infty}e^{i 2 \pi n t} \Phi(n)$$

and consider the case where $\phi(t) = e^{-a^2 t^2}$, $a \gt 0$. Then letting $u= a t$, we have

$$\sum_{n=-\infty}^{\infty} e^{-(u+n a)^2} = \frac{\sqrt{\pi}}{a} \left [1+2 \sum_{n=1}^{\infty} e^{-n^2 \pi^2/a^2} \cos{\left (2 \pi n \frac{u}{a} \right )} \right ]$$

The key observation here is that we can choose any value of $a$ we wish and this equation holds true. In fact, we can choose a value of $a$ such that the sum on the RHS may be ignored. Let's call this sum $\epsilon(u)$:

$$|\epsilon(u)| = 2 \left |\sum_{n=1}^{\infty} e^{-n^2 \pi^2/a^2} \cos{\left (2 \pi n \frac{u}{a} \right )}\right | \le \sum_{n=1}^{\infty} e^{-n^2 \pi^2/a^2}$$

Note that, when $a=1/2$ (which is the value used in the A & S formula above), $\left |\epsilon(u)\right | \le e^{-4 \pi^2} + e^{-16 \pi^2} + \cdots \approx 5 \cdot 10^{-17}$. Thus, we may rewrite the Poisson sum formula result as follows:

$$e^{u^2} [1+\epsilon(u)] = \frac{a}{\sqrt{\pi}} \left [1+2 \sum_{n=1}^{\infty} e^{-n^2 a^2} \cosh{2 n a u} \right ]$$

Now substitute this result into the integrals defining the error function above:

\begin{align}\frac{2}{\sqrt{\pi}}\int_0^y du \, e^{u^2} \sin{2 x u} &\approx \frac{2}{\sqrt{\pi}} \int_0^y du \frac{a}{\sqrt{\pi}} \left [1+2 \sum_{n=1}^{\infty} e^{-a^2 n^2} \cosh{2 n a u} \right ] \sin{2 x u}\\ &= 2 a \frac{1-\cos{2 x y}}{2 \pi x} + \frac{4 a}{\pi} \sum_{n=1}^{\infty} e^{-a^2 n^2} \int_0^y du \, \cosh{2 a n u}\, \sin{2 x u}\\ &= 2 a \frac{1-\cos{2 x y}}{2 \pi x} \\ &+ \frac{2 a}{\pi}\sum_{n=1}^{\infty} e^{-a^2 n^2} \frac{x (1-\cos{2 x y} \cosh{2 a n y})+ a n \sin{2 x y} \sinh{2 a n y}}{x^2+a^2 n^2} \end{align}

Similarly,

$$\frac{2}{\sqrt{\pi}}\int_0^y du \, e^{u^2} \cos{2 x u} \approx \\2 a \frac{\sin{2 x y}}{2 \pi x} + \frac{2 a}{\pi}\sum_{n=1}^{\infty} e^{-a^2 n^2} \frac{x \sin{2 x y} \cosh{2 a n y}+ a n \cos{2 x y} \sinh{2 a n y}}{x^2+a^2 n^2}$$

Putting this altogether, we reproduce the A & S result when $a=1/2$.

I have implemented this in Mathematica. One should note that the number of terms needed to reach a tolerance depends on the value of $z$, and is fairly sensitive to $\Im{z}$. Thus, I have implemented a simple while loop to achieve a desired precision.

It should be noted that the ceiling on this precision is the $10^{-16}$ rough figure I derived above. This, however, is of little importance, as this is the limit of what double precision, floating-point computation provides. That's why this result is a big deal: analytically, it is not equal to the error function, but computationally, it is equal for all practical purposes.

Anyway, here's the code:

 f[x_, y_, a_, n_] :=
Erf[x] + 2 a Exp[-x^2]/(2 Pi x) ((1 - Cos[2 x y]) +
I Sin[2 x y]) + (2 a Exp[-x^2]/Pi) Sum[
Exp[-a^2 k^2]/(x^2 +
a^2 k^2) ((x (1 - Cos[2 x y] Cosh[2 k a y]) +
a k Sin[2 x y] Sinh[2 a k y]) +
I (x Sin[2 x y] Cosh[2 k a y] +
a k Cos[2 x y] Sinh[2 a k y])), {k, 1, n}]


(This is the formula we derived above, in Mathematica syntax.)

g[x_, y_, a_] := Module[{n, err, tol, f1, f2, maxIters}, (
n = 1;
err = 1;
tol = 10^(-15);
maxIters = 50;
f2 = f[x, y, a, n];
While[Abs[err] > tol && n < maxIters, (
f1 = f2;
f2 = f[x, y, a, ++n];
err = Abs[f2/(f1 + 10^(-16)) - 1]
)];
f2
)] ;


(This is the while loop that computes erf to a desired precision. Note that the maxIters condition is necessary because there are points that seem to resist convergence. I think these may be zeroes of the error function, but I have not yet investigated.)

And now, here's a plot of some results; note that the plot of the effective number of decimal places of precision achieved.

ContourPlot[-Log[10,
Abs[g[x, y, 0.5]/(Erf[x + I y] + 10^(-16)) - 1]], {x, -2,
2}, {y, -4, 4}, PlotPoints -> 20, PlotLegends -> Automatic]


The high amount of detail is indicative of the noise we see when we have reached the limit of precision we can achieve. There is also some structure around where the computation was not able to achieve the desired level of precision; again, this is worth investigating.

• There is no $a$ on the LHS of your last approximation. – Pedro Tamaroff Mar 14 '14 at 21:29
• @PedroTamaroff: that's OK. Read through the derivation. – Ron Gordon Mar 14 '14 at 21:30
• Oh, my bad. =) ${}{}$ – Pedro Tamaroff Mar 14 '14 at 21:30
• I am bookmarking your answer! thanks, very helpful! I am trying to find some approximate solution to the integral from -inf to inf [f(t)*exp(i*(g(t) + cerf(t+id))] dt, where f(t), g(t) are some known functions, c and d are constants. I was thinking using the stationary phase method, and for that i need to separate the erf into real and imaginary. Hopefully using your comment and comment from GEdgar, I can get to something – Sleepyhead Mar 14 '14 at 21:34
• @Sleepyhead: wow, OK good luck. Note that you can use other values of $a$ than $1/2$ so you can tune the number of terms you need from the infinite sum. In that case, though, you need to re-estimate the max relative error. – Ron Gordon Mar 14 '14 at 22:04

Well, $$\text{Re}\;\text{erf}(a+ib) = \frac{\text{erf}(a+ib)+\text{erf}(a-ib)}{2}$$ and something similar for the imaginary part.

• This might work, thanks – Sleepyhead Mar 14 '14 at 20:33

The identities for the real and imaginary parts of the error function are:

\begin{aligned} {\mathop{\rm Re}\nolimits} \left[ {{\rm{erf}}\left( {x + iy} \right)} \right] &= \frac{{{\rm{erf}}\left( {x + iy} \right) + {\rm{erf}}\left( {x - iy} \right)}}{2}\\ {\mathop{\rm Im}\nolimits} \left[ {{\rm{erf}}\left( {x + iy} \right)} \right] &= \frac{{{\rm{erf}}\left( {x + iy} \right) - {\rm{erf}}\left( {x - iy} \right)}}{{2i}} \end{aligned}

See the paper [Abrarov and Quine, J. Math. Research, 7(2), 2015, pp. 163-174] for the rigorous proof: http://dx.doi.org/10.5539/jmr.v7n2