simpson's rule in C

simpsons.c

long double adaptive_simpsons(long double (*f)(long double), long double a, long double b, long double epsilon) {
	long double c = (a + b) / 2.0L;
	long double h = b - a;
	long double fa = f(a);
	long double fb = f(b);
	long double fc = f(c);
	long double integral = (h / 6.0L) * (fa + 4.0L * fc + fb);
	
	long double d = (a + c) / 2.0L;
	long double e = (c + b) / 2.0L;
	long double fd = f(d);
	long double fe = f(e);
	long double integral_left = (h / 12.0L) * (fa + 4.0L * fd + fc);
	long double integral_right = (h / 12.0L) * (fc + 4.0L * fe + fb);
	
	if( fabs(integral - integral_left - integral_right) < epsilon ) {
		return integral_left + integral_right;
	} else {
		return adaptive_simpsons(f, a, c, epsilon / 2.0L) + adaptive_simpsons(f, c, b, epsilon / 2.0L);
	}
}