數值積分 (辛普森公式 辛普森自适應法則) UVA

double F(double x) {
}
double simpson(double l, double r) {
	double m = (l + r) * 0.5;
	return (F(l) + 4 * F(m) + F(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double A) {
	double m = (l + r) * 0.5;
	double L = simpson(l, m), R = simpson(m, r);
	if (fabs(L + R - A) <= 15 * eps)
		return L + R + (L + R - A) / 15.0;
	return asr(l, m, eps * 0.5, L) + asr(m, r, eps * 0.5, R);
}
double asr(double l, double r, double eps) {
	return 2.0 * asr(l, r, eps, simpson(l, r));
}
int dblcmp(double x) {
	if (x < -eps)
		return -1;
	return x > eps;
}
           

UVA 1356

#include<cstdio>
#include<cstring>
#include<queue>
#include<vector>
#include<iostream>
#include<string>
#include<sstream>
#include<cctype>
#include<set>
#include<stack>
#include<numeric>
#include<functional>
#include<memory>
#include<deque>
#include<list>
#include<cmath>
#include<fstream>
#include<cstdlib>
#include<climits>
#include<iomanip>
#include<cstring>
#include<memory>
#include<bitset>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<VVI> VVVI;
const double eps = 1e-10;
const int maxn = 1111;
const int mod = 1000000007;
const int N = 1111;
const int M = 111111;
const int INF = 0x7f7f7f7f;
double a, w, h;
double F(double x) {
	return sqrt(1 + 4.0 * a * a * x * x);
}
double simpson(double l, double r) {
	double m = (l + r) * 0.5;
	return (F(l) + 4 * F(m) + F(r)) * (r - l) / 6.0;
}
double asr(double l, double r, double eps, double A) {
	double m = (l + r) * 0.5;
	double L = simpson(l, m), R = simpson(m, r);
	if (fabs(L + R - A) <= 15 * eps)
		return L + R + (L + R - A) / 15.0;
	return asr(l, m, eps * 0.5, L) + asr(m, r, eps * 0.5, R);
}
double asr(double l, double r, double eps) {
	return 2.0 * asr(l, r, eps, simpson(l, r));
}
int dblcmp(double x) {
	if (x < -eps)
		return -1;
	return x > eps;
}
int main() {
	double D, H, B, L;
	int cas;
	int n, T;
	double L1;
	scanf("%d", &T);
	for (cas = 1; cas <= T; ++cas) {
		scanf("%lf%lf%lf%lf", &D, &H, &B, &L);
		if (cas > 1)
			puts("");
		n = (int) ceil(B * 1.0 / D);
		D = B * 1.0 / n;
		L1 = L * 1.0 / n;
		double len, l, r, m;
		l = 0, r = H;
		w = D;
		while (r - l > eps) {
			m = (l + r) * 0.5;
			h = m;
			a = 4.0 * h / (w * w);
			len = asr(0, w * 0.5, 1e-6);
			if (dblcmp(len - L1) < 0)
				l = m;
			else
				r = m;
		}
		printf("Case %d:\n%.2lf\n", cas, 0.0 + H - h);
	}
	return 0;
}