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


double integrand(double t, void* params){
	double x = *(double*)params;
	return pow(t,x-1)*exp(-t);
}

double mygamma(double x){
	gsl_function F;
	F.function = &integrand;
	F.params = &x;
	const int limit = 1000;

	gsl_integration_workspace * w 
		= gsl_integration_workspace_alloc (limit);

	double r,e;
	int status = gsl_integration_qagiu(&F,0,1e-6,1e-6,limit,w,&r,&e);
	if(status!=0)fprintf(stderr,"mygamma: status=%i\n",status);

	gsl_integration_workspace_free (w);
	return r;
}