using System;
using static System.Math;
public static partial class math{

public delegate double integrand (double x);

public static double integrate
(integrand f, double a, double b, double acc, double eps){
	double f2=f(a+(b-a)/3), f3=f(a+2*(b-a)/3);
	double[] result=adapt(f,a,b,acc,eps,f2,f3);
	return result[0];
   }

private static double[] adapt
(integrand f,double a,double b,double acc,double eps,double f2,double f3){
	double h=b-a, f1=f(a+h/6), f4=f(a+5*h/6);
	double Q=(f1+f4)*h/3+(f2+f3)*h/6;
	double q=(f1+f2+f3+f4)*h/4;
	double tol=acc+eps*Abs(Q);
	double err=Abs(Q-q)/2;
	double[] result = {Q,err};
	if(err<Sqrt(2)*tol){return result;}
	else{
		double[] r1 = adapt(f,a,(a+b)/2,acc/Sqrt(2),eps,f1,f2);
		double[] r2 = adapt(f,(a+b)/2,b,acc/Sqrt(2),eps,f3,f4);
		double i=r1[0]+r2[0];
		err=Sqrt(r1[1]*r1[1]+r2[1]*r2[1]);
		double[] r = {i,err};
		return r;
      }
}//adapt

}//math