#include <stdio.h>
#include <math.h>
// GSL header:
#include "gsl/gsl_integration.h"

int ncalls=0;

// prototype of the adapt function:
double adapt(
	double (*f)(double),double a,double b,double eps,double acc,double *err);

// a function to integrate:
double f (double x, void * params) {
  double alpha = *(double*)params;
  double f = log(alpha*x) / sqrt(x);
  return f;
}

// a function to integrate:
double f2(double x){ncalls++; return log(x)/sqrt(x);}

int main (void) {

  double a=0; double b=1; double epsabs=0.0001; double epsrel=0.0001;
  size_t limit = 1000; double err,exact=-4;

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

  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);
  printf("acc,eps:          =  %g, %g\n",epsabs,epsrel);
  printf("exact             =  %.12f\n", exact);
  printf("---------------\nQAGS:\n");
  printf("result            =  %.12f\n", result);
  printf("estimated error   =  %.12f\n", error);
  printf("actual error      =  %.12f\n", fabs(exact-result));
  printf("ncalls            =  %d\n", w->size*21);

  epsrel=epsabs=1.e-8;
  double 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", fabs(exact-integ));
  printf("ncalls            = % d\n", ncalls);

  return 0;
}