#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include "debug.h"

double foo(double t, void* p){
	double result = exp(-t*t);
TRACE("foo: r=%g\n",result);
	return result;
}

double errorfc(double x){
	if(x<0)return 2-errorfc(-x);
	int limit=1000;
	gsl_function F = {&foo,NULL};
	double r,e;
	gsl_integration_workspace *w
		= gsl_integration_workspace_alloc (limit);
	gsl_integration_qagiu(&F,x,1e-6,1e-6,limit,w,&r,&e);
	return r*2/sqrt(M_PI);
	gsl_integration_workspace_free (w); // !!!
}