#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 *);

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);}
double f2(double x){ncalls++; return log(x)/sqrt(x);}

int main (void) {

  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  double result, error; double alpha = 1.0;

double a=0; double b=1; double epsabs=0.001; double epsrel=0.001;
size_t limit=1000; int key=1; double integ,err,exact;

  gsl_function F; F.function = &f; F.params = &alpha;
  gsl_integration_qags(&F,a,b,epsabs,epsrel,limit,w,&result,&error); 

//  exact=atan(b)-atan(a);
  exact=-4;
  printf("acc,eps:          =  %g, %g\n",epsabs,epsrel);
  printf("exact             =  %.12f\n", exact);
  printf("---------------\nQAG:\n");
  printf("result            =  %.12f\n", result);
  printf("estimated error   =  %.12f\n", error);
  printf("actual error      =  %.12f\n", exact-result);
  printf("ncalls            =  %d\n", w->size*21);

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

  printf("---------------\nadapt24:\n");
  printf("result            = % .12f\n", integ);
  printf("estimated error   = % .12f\n", err);
  printf("actual error      =  %.12f\n", exact-integ);
  printf("ncalls            = % d\n", ncalls);

  return 0;
}