#define real double
#include <math.h>

real adapt24(real (*func)(real),
	real a,real b,real f2, real f3, real eps,real acc,real *err){

	real h=b-a; real f1=(*func)(a+h/6); real f4=(*func)(a+5*h/6);
	real j4=(f1+f4)*h/3+(f2+f3)*h/6; real j2=(f1+f2+f3+f4)*h/4;
	real tol=acc+eps*fabs(j4); *err=fabs(j4-j2)/7;

	if(*err < sqrt(2.0)*tol) return j4;
	else	{
		real err1,err2;
		real i1=adapt24(func,a,(a+b)/2,f1,f2,eps,acc/sqrt(2.),&err1);
		real i2=adapt24(func,(a+b)/2,b,f3,f4,eps,acc/sqrt(2.),&err2);
		*err=sqrt(err1*err1+err2*err2);
		return i1+i2; } }

real adapt(real (*func)(real),real a,real b,real eps,real acc,real *err){
	real f2=(*func)(a+(b-a)/3); real f3=(*func)(a+2*(b-a)/3);
	real integ=adapt24(func,a,b,f2,f3,eps,acc,err);
	return integ; }