Least Action

Nontrivializing triviality..and vice versa.

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

with 5 comments

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,

erf(z) = \frac{2}{\sqrt{\pi}}\sum_{n = 0}^{\infty}\frac{(-1)^n z^{2n+1}}{n!(2n+1)}

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;
}
Advertisements

Written by Vivek

January 28, 2011 at 21:12

5 Responses

Subscribe to comments with RSS.

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

    stevengj

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

      Vivek

      November 11, 2013 at 01:38

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

    joachimwuttke

    May 20, 2013 at 17:21

    • Thanks for the reference!

      Vivek

      November 11, 2013 at 01:35

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

    Joachim Wuttke

    November 2, 2013 at 19:41


Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: