chyx's Blog

http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method

double simpson_rec( double a, double b, double S, double fa, double fb, double fc,
		int bottom, double epsilon){
	double c = (a+b)/2, h = b-a;
	double d = (a+c)/2, e = (c+b)/2;
	double fd = f(d), fe = f(e);
	double SL = (h/12)*(fa + 4*fd+fc);
	double SR = (h/12)*(fc + 4*fe+fb);
	double S2 = SL + SR;
	if( bottom <= 0 || fabs(S2 - S ) <= 15*epsilon){
		return S2 + (S2-S)/15;
	}
	return 	simpson_rec(a,c,SL,fa,fc,fd, bottom-1, epsilon/2) + 
		simpson_rec(c,b,SR,fc,fb,fe, bottom-1, epsilon/2);
}

double simpson( double a, double b, double epsilon, int rec){
	double c = (a+b)/2, h = b-a;
	double fa = f(a);
	double fb = f(b);
	double fc = f(c);
	double S = (h/6) * (fa+4*fc+fb);
	return simpson_rec(a,b,S,fa,fb,fc,rec,epsilon);
}