# Least Action

Nontrivializing triviality..and vice versa.

## 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,

$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);
}

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


Written by Vivek

January 28, 2011 at 21:12

## CUDA on Ubuntu Maverick Meerkat 10.10

with one comment

(To be edited)

If you’re getting errors like

/usr/bin/ld: cannot find -lGL
/usr/bin/ld: cannot find -lGLU
/usr/bin/ld: cannot find -lX11
/usr/bin/ld: cannot find -lXi
/usr/bin/ld: cannot find -lXmu
/usr/bin/ld: cannot find -lglut


then you need to do the following

sudo apt-get install libxi libxi-dev

Running make again will result in the following errors:

/usr/bin/ld: cannot find -lGL
/usr/bin/ld: cannot find -lGLU
/usr/bin/ld: cannot find -lXmu
/usr/bin/ld: cannot find -lglut


Now, let’s execute:

sudo apt-get install freeglut3 freeglut3-dev

This brings down the errors to

/usr/bin/ld: cannot find -lXmu

So, we just have to do one more apt-get:

sudo apt-get install libxmu6 libxmu-dev

—-
If you get errors suggesting that your libcudart.so.1 is missing, it means your LD_CONFIG_PATH isn’t set right. To set it permanently, use

sudo ldconfig -v /usr/local/cuda/lib64/

on a 64-bit system and

sudo ldconfig -v /usr/local/cuda/lib/

on a 32 bit system.

Written by Vivek

January 25, 2011 at 22:49