Fast Complex Error Function and Faddeeva Function

While the real error function is easy to approximate and is available in C++ standard library, the imaginary and complex valued error functions are neither.

Fig 1. Logarithmic plot of the absolute value of error function in the complex plane.

The complex error function as well as the complementary and imaginary error functions

$$ \begin{aligned} \textrm{erf}\left(z\right) &= \frac{2}{\pi} \int_0^z e^{-\xi^2}d\xi \\ \textrm{erfc}\left(z\right) &= 1-\textrm{erf}\left(z\right) \\ \textrm{erfi}\left(z\right) &= -i\textrm{erf}\left(iz\right) \end{aligned} $$

are related to the Faddeeva function $w(z)$: $$ w\left(z\right) = e^{-z^2}\textrm{erfc}\left(-iz\right) $$ and therefore $$ \textrm{erf}\left(z\right) = 1 - e^{-z^2}w\left(iz\right) $$

The Faddeeva function can be approximated in the upper half of the complex plane ($\Im\left\{z\right\}\ge 0$) with configurable precision using an iterative expansion [1]. The constant iteration count makes it very suitable for GPU work and by using $N=8$ precomputed terms a very fast implementation is realized. For example, in GLSL, using single-precision complex arithmetics:

// The Faddeeva function
vec2 faddeeva(vec2 z) {
	const int M = 4;
	const int N = 1 << (M-1);

	// Precomputed Faddeeva function approximation coefficients (M==4)
	const vec2 A[N] = {
		vec2(+0.983046454995208, 0.0),
		vec2(-0.095450491368505, 0.0),
		vec2(-0.106397537035019, 0.0),
		vec2(+0.004553979597404, 0.0),
		vec2(-0.000012773721299, 0.0),
		vec2(-0.000000071458742, 0.0),
		vec2(+0.000000000080803, 0.0),
		vec2(-0.000000000000007, 0.0),
	};
	const vec2 B[N] = {
		vec2(0.0, -1.338045597353875),
		vec2(0.0, +0.822618936152688),
		vec2(0.0, -0.044470795125534),
		vec2(0.0, -0.000502542048995),
		vec2(0.0, +0.000011914499129),
		vec2(0.0, -0.000000020157171),
		vec2(0.0, -0.000000000001558),
		vec2(0.0, +0.000000000000003),
	};
	const vec2 C[N] = {
		vec2(0.392699081698724, 0.0),
		vec2(1.178097245096172, 0.0),
		vec2(1.963495408493621, 0.0),
		vec2(2.748893571891069, 0.0),
		vec2(3.534291735288517, 0.0),
		vec2(4.319689898685965, 0.0),
		vec2(5.105088062083414, 0.0),
		vec2(5.890486225480862, 0.0),
	};
	const float s = 2.75;

	// Constrain to imag(z)>=0
	const float sgni = imag(z)<0 ? -1 : 1;
	z *= sgni;

	// Approximate
	const vec2 t = z + vec2(0,0.5)*s;
	vec2 w = vec2(.0);
	for (int m=0; m<N; ++m)
		w += cdiv(A[m] + cprod(t,B[m]), csqr(C[m]) - csqr(t));

	// Invert back
	if (sgni < 0)
		w = 2.0*cexp(-csqr(z)) - w;

	return w;
}
// Error function
vec2 erf(vec2 z) {
	vec2 z_1i = cprod(vec2(0,1), z);	// 1i*z
	return vec2(1,0) - cprod(cexp(-csqr(z)), faddeeva(z_1i));
}

With the complex arithmetic functions:

float real(vec2 z) 				{ return z.x; }
float imag(vec2 z) 				{ return z.y; }
vec2  cprod(vec2 a, vec2 b) 	{ return vec2(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
vec2  csqr(vec2 a) 				{ return vec2(a.x*a.x-a.y*a.y, 2*a.x*a.y); }
vec2  cdiv(vec2 a, vec2 b) 	{ return vec2(a.x*b.x+a.y*b.y, a.y*b.x-a.x*b.y) * (1.0/(b.x*b.x+b.y*b.y)); }
vec2  cexpi(float d) 			{ return vec2(cos(d), sin(d)); }
vec2  cexp(vec2 z) 				{ return exp(z.x) * cexpi(z.y); }

References

[1] S.M. Abrarov, B.M. Quine, Efficient algorithmic implementation of the Voigt/complex error function based on exponential series approximation, Applied Mathematics and Computation, Volume 218, Issue 5, 2011, Pages 1894-1902, ISSN 0096-3003, https://doi.org/10.1016/j.amc.2011.06.072