#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#define real double

int ncalls=0;

real adapt(real (*)(real),real,real,real,real,real,real,real *);

double f (double x, void * params) {
  double alpha = *(double *) params;
//  double f = log(alpha*x) / sqrt(x);
  double f=1/(1+x*x);
  return f;
}

double f2(double x){ncalls++; return 1/(1+x*x);}

int main (void)
{
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  
  double result, error; double alpha = 1.0;

double a=-20; double b=20; double epsabs=0.0001; double epsrel=0.0001;
size_t limit=1000; int key=1; double integ,err,fa,fb;

  gsl_function F;
  F.function = &f;
  F.params = &alpha;

  gsl_integration_qags(&F,a,b,epsabs,epsrel,limit,w,&result,&error); 

  printf("acc,eps:          =  %g, %g\n",epsabs,epsrel);
  printf("QAG:\n");
  printf("result            =  %.12f\n", result);
  printf("estimated error   =  %.12f\n", error);
  printf("ncalls            =  %d\n", w->size*21);

  fa=f2(a); fb=f2(b);
  integ=adapt(&f2,a,b,fa,fb,epsrel,epsabs,&err);

  printf("adapt:\n");
  printf("result            = % .12f\n", integ);
  printf("estimated acc,eps = % .12f, %.12f\n", err,err/integ);
  printf("ncalls            = % d\n", ncalls);

  return 0;
}