[BZOJ 2178] 圆的面积并 【Simpson积分】

题目链接：BZOJ - 2178

题目分析

用Simpson积分，将圆按照 x 坐标分成连续的一些段，分别用 Simpson 求。

注意：1)Eps要设成 1e-13  2）要去掉被其他圆包含的圆。

代码

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;

typedef double LF;

const LF Eps = 1e-13;

const int MaxN = 1000 + 5;

int n, Lc, Rc, Top, Tot;

LF Lx, Rx, Ans;

inline LF gmin(LF a, LF b) {return a < b ? a : b;}
inline LF gmax(LF a, LF b) {return a > b ? a : b;}
inline LF Sqr(LF x) {return x * x;}

struct Point
{
	LF x, y;
	Point() {}
	Point(LF a, LF b)
	{
		x = a; y = b;
	}
};

inline LF Dis(Point p1, Point p2)
{
	return sqrt(Sqr(p1.x - p2.x) + Sqr(p1.y - p2.y));
}

struct Circle
{
	Point o;
	LF r;
} C[MaxN], S[MaxN];

inline bool Cmp1(Circle c1, Circle c2)
{
	return c1.r < c2.r;
}

inline bool Cmp2(Circle c1, Circle c2)
{
	return c1.o.x - c1.r < c2.o.x - c2.r;
}

bool Del[MaxN];

struct Segment
{
	LF l, r;
} Seg[MaxN];

inline bool Cmp3(Segment s1, Segment s2)
{
	return s1.l < s2.l;
}

inline LF f(LF x)
{
	LF ret = 0.0, p, q, Hi;
	Top = 0;
	for (int i = Lc; i <= Rc; ++i)
	{
		if (x <= S[i].o.x - S[i].r || x >= S[i].o.x + S[i].r) continue;
		Hi = sqrt(Sqr(S[i].r) - Sqr(S[i].o.x - x));
		Seg[++Top].l = S[i].o.y - Hi; Seg[Top].r = S[i].o.y + Hi;
	}
	sort(Seg + 1, Seg + Top + 1, Cmp3);
	for (int i = 1, j; i <= Top; ++i)
	{
		p = Seg[i].l; q = Seg[i].r;
		for (j = i + 1; j <= Top; ++j)
		{
			if (Seg[j].l > q) break;
			if (Seg[j].r > q) q = Seg[j].r;
		}
		ret += q - p; i = j - 1;	
	}
	return ret;	
}

inline LF Simpson(LF l, LF r, LF fl, LF fmid, LF fr)
{
	return (fl + 4.0 * fmid + fr) * (r - l) / 6.0;
}

inline LF RSimpson(LF l, LF r, LF fl, LF fmid, LF fr)
{
	LF mid, p, q, x, y, z;
	mid = (l + r) / 2.0;
	p = f((l + mid) / 2.0); q = f((mid + r) / 2.0);
	x = Simpson(l, r, fl, fmid, fr);
	y = Simpson(l, mid, fl, p, fmid);
	z = Simpson(mid, r, fmid, q, fr);
	if (fabs(x - y - z) < Eps) return y + z;
	else return RSimpson(l, mid, fl, p, fmid) + RSimpson(mid, r, fmid, q, fr);
}

int main()
{
	scanf("%d", &n);
	int a, b, c;
	for (int i = 1; i <= n; ++i)
	{
		scanf("%d%d%d", &a, &b, &c);		
		C[i].o = Point((LF)a, (LF)b);
		C[i].r = (LF)c;
	}	
	sort(C + 1, C + n + 1, Cmp1);
	memset(Del, 0, sizeof(Del));
	for (int i = 1; i <= n; ++i)
		for (int j = i + 1; j <= n; ++j)
			if (Dis(C[i].o, C[j].o) <= C[j].r - C[i].r)
			{
				Del[i] = true;
				break;
			}
	Tot = 0;
	for (int i = 1; i <= n; ++i) 
		if (!Del[i]) S[++Tot] = C[i];
	sort(S + 1, S + Tot + 1, Cmp2);
	Ans = 0.0;
	for (int i = 1, j; i <= Tot; ++i)
	{
		Lc = i; Rc = i; Lx = S[i].o.x - S[i].r; Rx = S[i].o.x + S[i].r;
		for (j = i + 1; j <= Tot; ++j)
		{
			if (S[j].o.x - S[j].r > Rx) break;
			if (S[j].o.x + S[j].r > Rx) Rx = S[j].o.x + S[j].r;
		}
		Rc = j - 1; i = j - 1;
		Ans += RSimpson(Lx, Rx, f(Lx), f((Lx + Rx) / 2.0), f(Rx));
	}
	printf("%.3lf
", Ans);
	return 0;
}