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

int n=0;
double integrand (double x, void* params){
	n++;
	return sqrt(1-x*x);
}

double mypi(void){
	int limit=1000;
//	gsl_integration_workspace* w = gsl_integration_workspace_alloc(limit);
	gsl_integration_cquad_workspace* w = gsl_integration_cquad_workspace_alloc(limit);
	gsl_function f;
	f.function = integrand;
	f.params=NULL;

	double result, err;
//	int status=gsl_integration_qags (&f, 0, 1, 1e-12, 0, limit, w, &result, &err);
	size_t nevals;
	int status=gsl_integration_cquad (&f, 0, 1, 1e-12, 0, w, &result, &err, &nevals);

	if(status !=0) fprintf(stderr,"mypi: status=%i\n",status);
	fprintf(stderr,"mypi: nevals = %i\n",(int)nevals);
	fprintf(stderr,"mypi: err = %g\n",err);

	gsl_integration_cquad_workspace_free(w);

	return 4*result;
}

int main(){
	double pi = mypi();
	printf("pi=\t%20.15g\n",pi);
	printf("PI=\t%20.15g\n",M_PI);
	printf("ratio=\t%20.15g\n",pi/M_PI);
	printf("n=%i\n",n);
return 0;
}