## Evaluation of Complex Error Functions erf(z) using GSL/C

For a small QFT calculation, I needed to numerically evaluate the imaginary error function erfi(x) = erf(i x). But it turns out that GSL (and most other numerical recipe code I could find) can only deal with erf(x), where x is real. Here’s a poor man’s implementation of erf(z) through a standard Taylor expansion,

The catch here is to deal with propagation of errors in the *complex* Taylor series, and to also somehow benchmark the results. Well, I haven’t been able to think about this yet, but I was able to confirm that for **real** arguments, my Taylor series code is about as good as the GSL error function gsl_sf_erf(x).

I am working on a CUDA implementation of this now, because in my project, I need to perform a numerical integration over the error function, which is quite intensive even for Mathematica. For now, I’m just sharing the serial implementation. (**Hint**: if you can write wrapper for each of the gsl functions inside my Taylor series calculating method, you can embed them within a __global__ kernel call through a struct. Alternatively — and less fun — you can just parallelize the for loop.)

/* Function to compute erf(z) using a Taylor series expansion /* Author: Vivek Saxena /* Last updated: January 28, 2011 21:11 hrs */ #include <gsl/gsl_complex_math.h> #include <gsl/gsl_sf_erf.h> #include <gsl/gsl_sf_pow_int.h> #include <gsl/gsl_sf_gamma.h> #include <stdio.h> #define PI 3.1415926543 double cz = 2/sqrt(PI); const int TERMS = 10; // no of terms to use in the Taylor series gsl_complex erfTaylor( gsl_complex z, int trunc ){ gsl_complex res = gsl_complex_rect(0,0), num = gsl_complex_rect(0,0), den = gsl_complex_rect(1,0), snum = gsl_complex_rect(1,0), temp = gsl_complex_rect(0,0); for(int i = 0; i < trunc; i++){ snum = gsl_complex_rect( cz * gsl_sf_pow_int(-1, i), 0 ); num = gsl_complex_mul(snum, gsl_complex_pow_real(z, 2*i+1)); den = gsl_complex_rect((2*i + 1)*gsl_sf_fact(i),0); temp = gsl_complex_div(num, den); res = gsl_complex_add(res, temp); } return res; } int main ( void ){ printf( "Real error function\n\n"); for ( float i = 0; i <= 1; i += 0.01 ){ float gslerror = gsl_sf_erf(i); float taylor = GSL_REAL( erfTaylor( gsl_complex_rect(i, 0), 10 ) ); printf("erf(%f): gsl = %f, taylor = %f, mag error = %f\n", i, gslerror, taylor, abs(gslerror-taylor)); } gsl_complex t, arg; printf( "\n\nImaginary error function\n\n"); for (float i = 0; i <= 1; i += 0.01 ){ /* this would be your generic argument z in erf(z). * so if z = x + iy, then * arg = gsl_complex_rect(x, y); */ arg = gsl_complex_rect(0, i); t = erfTaylor( arg, TERMS ); printf("erf(%f + i %f) = %f + i %f\n", GSL_REAL(arg), GSL_IMAG(arg), GSL_REAL(t), GSL_IMAG(t)); } return 0; }

Unfortunately, the Taylor series is a very inefficient way to compute erf(z) when |z| is not small. (Note also that you have some gratuitous inaccuracy because you give pi to only 11 digits.) For example, with z=1+2i and using 10 terms as above, your code gives -0.273745-4.74766i when the correct answer is -0.536644…-5.04914…i. With 20 terms, your code gives the correct answer, but to only 6 digits.

Fortunately, there is an alternative: I have a free/open-source C/C++ routine that computes erf(z) to nearly machine precision for all complex z, available at: http://ab-initio.mit.edu/Faddeeva … besides being more accurate, on my machine it is about 5 times faster for z=1+2i than your code above with 10 terms.

stevengjDecember 23, 2012 at 10:42

Thanks Steven, I think I used your notes on Perfectly Matched Layers about 3 years ago for my Bachelor’s thesis. I thoroughly benefited from them, and I am grateful to you for bringing this to my notice. At the time I wrote this, this was a quick and dirty way of getting some work done which did not depend on a very good value of erf(z). And I conveniently forgot about it, but your comment will prove to be useful to those who may stumble upon this post and rely on its numerical accuracy.

VivekNovember 11, 2013 at 01:38

for a mathematically sound and userfriendly implementation of erf(z), see code written by Steven G. Johnson and packaged by myself: http://apps.jcns.fz-juelich.de/doku/sc/libcerf. Regards – Joachim Wuttke

joachimwuttkeMay 20, 2013 at 17:21

Thanks for the reference!

VivekNovember 11, 2013 at 01:35

For a full-fledged numerical library providing complex error functions and related functions, see my package libcerf, apps.jcns.fz-juelich.de/libcerf.

Joachim WuttkeNovember 2, 2013 at 19:41