自適應辛普森法

自適應辛普森法

double a, b, c, d, l, r;

const double eps = 1e-8;
double F(double x) {
    //需要積分的公式
	return (c * x + d) / (a * x + b);
}
double simpson(double l, double r) {
	double mid = (l + r) / 2;
	return (F(l) + 4 * F(mid) + F(r)) * (r - l) / 6;
}
double asr(double l, double r, double A) {
	double mid = (l + r) / 2;
	double L = simpson(l, mid), R = simpson(mid, r);
	if (fabs(L + R - A) <= 15 * eps)
		return L + R + (L + R - A) / 15.0;
	return asr(l, mid, L) + asr(mid, r, R);
}
double asr(double l, double r) {
	return asr(l, r, simpson(l, r));
}