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.
The complex error function as well as the complementary and imaginary error functions
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
