Rockdu's Blog

【题目描述】

给出N个圆，求其面积并

【输入】

先给一个数字N ,N< = 1000

接下来是N行是圆的圆心，半径，其绝对值均为小于1000的整数

【输出】

面积并，保留三位小数

辛普森积分实在太强了……

它相当于用很多1、2、3次函数去拟合一个函数以达到轻松计算积分的目的。

这里给一篇博客，讲的很好：一篇博客

然后有一些小小的优化：
为了解决函数不连续的问题，可以设置一个限定长度，让辛普森积分只在值域范围在限定长度以内的情况下才停止递归。

然后有一个精度优化，辛普森的eps可以设置成每向下递归一层就除以2，效果拔群。

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;
const int maxn = 1e3 + 5;
struct Circle {
	double x, y, R;
}cir[maxn];
struct Itv {
	double b, t;
	bool operator < (const Itv & A) const {
		return b < A.b;
	}
}itv[maxn];
int n, cnt;
double dx, tmp, L = 1e300, R = -1e300;
#define sqr(x) ((x) * (x))
double F(double x) {
	cnt = 0;
	double gu, mxT = -1e300, mxB = -1e300, ans = 0;
	for(register int i = 1; i <= n; ++i) {
		if((tmp = abs(cir[i].x - x)) > cir[i].R) continue;
		gu = sqrt(sqr(cir[i].R) - sqr(tmp));
		itv[++cnt] = (Itv){cir[i].y - gu, gu + cir[i].y};
	}
	sort(itv + 1, itv + cnt + 1);
	for(register int i = 1; i <= cnt; ++i) {
		if(mxT < itv[i].b) ans += mxT - mxB, mxB = itv[i].b, mxT = itv[i].t;
		else mxT = max(mxT, itv[i].t);
	}
	ans += mxT - mxB;
	return ans;
}

double S(double L, double R, double FL, double FR, double FM) {
	return (R - L) * (FL + 4 * FM + FR) / 6.0;
}

double Simpson(double L, double R, double FL, double FM, double FR, double eps) {
	double mid = (L + R) / 2.0, nans = S(L, R, FL, FR, FM);
	double FLM = F((L + mid) / 2.0), FRM = F((mid + R) / 2.0);
	if(R - L <= 1)
		if(abs(S(L, mid, FL, FM, FLM) + S(mid, R, FM, FR, FRM) - nans) <= 15 * eps) 
			return nans;
	nans = Simpson(L, mid, FL, FLM, FM, eps / 2) + Simpson(mid, R, FM, FRM, FR, eps / 2);
	return nans;
}

int main() {
	scanf("%d", &n);
	for(register int i = 1; i <= n; ++i) {
		scanf("%lf%lf%lf", &cir[i].x, &cir[i].y, &cir[i].R);
		L = min(L, cir[i].x - cir[i].R);
		R = max(R, cir[i].x + cir[i].R);
	}
	printf("%.3lf", Simpson(L, R, F(L), F((L + R) / 2.0), F(R), 1e-5));
	return 0;
}