// (C) 2020 Dmitri Fedorov; License: GNU GPL v3+; no warranty.
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){
/// attempts to calculate int_a^b f(x)dx
/// with absolute precision acc or relative precision eps
	//return adapto4(f,a,b,acc,eps);
	//double f1=f(a), f3=f((a+b)/2), f5=f(b);
	//return adaptc5(f,a,b,acc,eps,f1,f3,f5);
	return adaptc7(f,a,b,acc,eps);
}//integrate

private static double adapto4
(integrand f,double a,double b,double acc,double eps,
double f2=double.NaN,double f3=double.NaN){
/// four point open adaptive integrator
	double h=b-a, f1=f(a+h/6), f4=f(a+5*h/6), sqr2=Sqrt(2);
	if(double.IsNaN(f2)){f2=f(a+2*h/6);f3=f(a+4*h/6);}
	double approx1=(3*f1+4*f2      +5*f4)*h/12;
	double approx2=(5*f1      +4*f3+3*f4)*h/12;
	double integral=(approx1+approx2)/2;
	double error=Abs(approx1-approx2)/2;
	double tolerance=acc+eps*Abs(integral);
	if(error<sqr2*tolerance)
		return integral;
	else
		return 	adapto4(f,a,(a+b)/2,acc/sqr2,eps,f1,f2)+
			adapto4(f,(a+b)/2,b,acc/sqr2,eps,f3,f4);
}//adapto

private static double adaptc5
(integrand f,double a,double b,double acc,double eps,double f1,double f3,double f5)
{ /// five point closed adaptive integrator
	double h=b-a, f2=f(a+h/4), f4=f(a+3*h/4), sqr2=Sqrt(2);
	double integral=(7*f1+32*f2+12*f3+32*f4+7*f5)*h/90;
	double approx  =(  f1+ 8*f2+       8*f4+  f5)*h/18;
	double error=Abs(integral-approx)/2;
	double tolerance=acc+eps*Abs(integral);
	if(error<sqr2*tolerance)
		return integral;
	else
		return 	adaptc5(f,a,(a+b)/2,acc/sqr2,eps,f1,f2,f3)+
			adaptc5(f,(a+b)/2,b,acc/sqr2,eps,f3,f4,f5);
}//adaptc5

private static double adaptc7
(integrand f,double a,double b,double acc,double eps,
double f1=double.NaN,double f3=double.NaN,double f5=double.NaN,double f7=double.NaN)
{ /// five point closed adaptive integrator
	double h=b-a,f2=f(a+h/6),f4=f(a+3*h/6),f6=f(a+5*h/6),sqr2=Sqrt(2);
	if(double.IsNaN(f1)){f1=f(a);f3=f(a+2*h/6);f5=f(a+4*h/6);f7=f(b);}
	double integral=
(41*f1/840+9*f2/35+9*f3/280+34*f4/105+9*f5/280+9*f6/35+41*f7/840)*h;
	double approx  =
(13*f1/200+4*f2/25+11*f3/40+          11*f5/40+4*f6/25+13*f7/200)*h;
	double error=Abs(integral-approx)/2;
	double tolerance=acc+eps*Abs(integral);
	if(error<sqr2*tolerance)
		return integral;
	else
		return 	adaptc7(f,a,(a+b)/2,acc/sqr2,eps,f1,f2,f3,f4)+
			adaptc7(f,(a+b)/2,b,acc/sqr2,eps,f4,f5,f6,f7);
}//adaptc5

}//math