自适应simpson积分，刚开始写了个朴素的，结果 WA on #9 ，后来又减小了一点误差，但是还是不对。非递归拟合也不对。我直接爬了。。有没有巨神帮忙瞧瞧？

• 朴素 simpson

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define rep(i, a, b) for (int i = (a); i <= (b); i ++ )
#define rop(i, a, b) for (int i = (a); i < (b); i ++ )
#define dep(i, a, b) for (int i = (a); i >= (b); i -- )
#define dop(i, a, b) for (int i = (a); i > (b); i -- )

using namespace std;

using LL = long long;
using PII = pair<int, int>;
using PLL = pair<LL, LL>;
using PDD = pair<double, double>;

const int N = 510;
const double eps = 1e-8;

struct S {
	double h, r0, r1;
	double x0, x1;
	double s0, s1;
	double k, b;
	void get() {
		if (fabs(r0 - r1) < eps) {
			s0 = x0, s1 = x1;
			k = 0, b = r0;
			return;
		}
		if (r0 > r1) {
//			cout << (double)(r0 - r1) / (x1 - x0) << endl;
			double theta = asin((double)(r0 - r1) / (x1 - x0));
//			cout << theta << endl;
			k = - tan(theta), s0 = x0 + r0 * sin(theta);
			s1 = x1 + r1 * sin(theta), b = r0 * cos(theta) - k * s0;
		}
		if (r0 < r1) {
			double theta = asin((double)(r0 - r1) / (x0 - x1));
			k = tan(theta), s0 = x0 - r0 * sin(theta);
			s1 = x1 - r1 * sin(theta), b = r0 * cos(theta) - k * s0;
		}
	}
}p[N];
int n;
double alpha;

double sqr(double x) { return x * x; }
double f(double x) {
	double maxn = 0;
	for (int i = 1; i <= n; i ++ ) {
		if (x >= p[i].s0 && x <= p[i].s1)
			maxn = max(maxn, p[i].k * x + p[i].b);
		if (x >= p[i].x0 - p[i].r0 && x <= p[i].x0 + p[i].r0)
			maxn = max(maxn, sqrt(sqr(p[i].r0) - sqr(p[i].x0 - x)));
	}
	return maxn;
}

double simpson(double l, double r) {
	double mid = (l + r) / 2.0;
	return (r - l) * (f(r) + f(l) + 4 * f(mid)) / 6;
}

double query(double l, double r, double area) {
	double mid = (l + r) / 2.0;
	double left = simpson(l, mid), right = simpson(mid, r);
	if (fabs(left + right - area) < eps) return left + right;
	else return query(l, mid, left) + query(mid, r, right);
}

int main() {
	scanf("%d%lf", &n, &alpha);
	double c = (double)1 / tan(alpha);
	scanf("%lf", &p[0].h); // the Whole Height Of the Tree
	p[0].x1 = p[0].h * c;
	for (int i = 1; i <= n; i ++ )
		scanf("%lf", &p[i].h);
	for (int i = 1; i <= n; i ++ ) {
		scanf("%lf", &p[i].r0);
//		p[i].r0 /= 2;
	}
	
	double l = 1e10, r = -1e10;
	for (int i = 1; i <= n; i ++ ) {
		p[i].x0 = p[i - 1].x1, p[i].x1 = p[i].x0 + p[i].h * c;
		p[i].r1 = p[i + 1].r0;
		p[i].get();
		l = min(l, p[i].x0 - p[i].r0);
		r = max(l, p[i].x0 + p[i].r0);
	}
	l = min(l, p[1].x0 - p[1].r0);
	r = max(r, p[n].x1);
	
	printf("%.2lf\n", query(l, r, simpson(l, r)) * 2.0);
	return 0;
}

• 高拟合 simpson

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define rep(i, a, b) for (int i = (a); i <= (b); i ++ )
#define rop(i, a, b) for (int i = (a); i < (b); i ++ )
#define dep(i, a, b) for (int i = (a); i >= (b); i -- )
#define dop(i, a, b) for (int i = (a); i > (b); i -- )

using namespace std;

using LL = long long;
using PII = pair<int, int>;
using PLL = pair<LL, LL>;
using PDD = pair<double, double>;

const int N = 510;

struct S {
	double h, r0, r1;
	double x0, x1;
	double s0, s1;
	double k, b;
	void get() {
		if (r0 == r1) {
			s0 = x0, s1 = x1;
			k = 0, b = r0;
			return;
		}
		if (r0 > r1) {
//			cout << (double)(r0 - r1) / (x1 - x0) << endl;
			double theta = asin((double)(r0 - r1) / (x1 - x0));
//			cout << theta << endl;
			k = - tan(theta), s0 = x0 + r0 * sin(theta);
			s1 = x1 + r1 * sin(theta), b = r0 * cos(theta) - k * s0;
		}
		if (r0 < r1) {
			double theta = asin((double)(r0 - r1) / (x0 - x1));
			k = tan(theta), s0 = x0 - r0 * sin(theta);
			s1 = x1 - r1 * sin(theta), b = r0 * cos(theta) - k * s0;
		}
	}
}p[N];
int n;
double alpha;

double sqr(double x) { return x * x; }
double f(double x) {
	double maxn = 0;
	for (int i = 1; i <= n; i ++ ) {
		if (x >= p[i].s0 && x <= p[i].s1)
			maxn = max(maxn, p[i].k * x + p[i].b);
		if (x >= p[i].x0 - p[i].r0 && x <= p[i].x0 + p[i].r0)
			maxn = max(maxn, sqrt(sqr(p[i].r0) - sqr(p[i].x0 - x)));
	}
	return maxn;
}

double simpson(double l, double r) {
	double mid = (l + r) / 2.0;
	return (r - l) * (f(r) + f(l) + 4 * f(mid)) / 6;
}

double query(double l, double r, double eps, double area, double step) {
	double mid = (l + r) / 2.0;
	double left = simpson(l, mid), right = simpson(mid, r);
	if (fabs(left + right - area) < eps * 15 && step < 0)
		return left + right + (left + right - area) / 15;
	return query(l, mid, eps / 2.0, left, step - 1) + query(mid, r, eps / 2.0, right, step - 1);
}

int main() {
	scanf("%d%lf", &n, &alpha);
	double c = (double)1 / tan(alpha);
	scanf("%lf", &p[0].h); // the Whole Height Of the Tree
	p[0].x1 = p[0].h * c;
	for (int i = 1; i <= n; i ++ )
		scanf("%lf", &p[i].h);
	for (int i = 1; i <= n; i ++ ) {
		scanf("%lf", &p[i].r0);
//		p[i].r0 /= 2;
	}
	
	double l = 1e10, r = -1e10;
	for (int i = 1; i <= n; i ++ ) {
		p[i].x0 = p[i - 1].x1, p[i].x1 = p[i].x0 + p[i].h * c;
		p[i].r1 = p[i + 1].r0;
		p[i].get();
		l = min(l, p[i].x0 - p[i].r0);
		r = max(l, p[i].x0 + p[i].r0);
	}
	l = min(l, p[1].x0 - p[1].r0);
	r = max(r, p[n].x1);
	
	printf("%.2lf\n", query(l, r, 0.0001, simpson(l, r), 15) * 2.0);
	return 0;
}