BZOJ2178 圆的面积并

计算几何题。别人都用simpson搞。然后我决定不水一发，用hwd神犇在wc上讲的关键点法做。

然后就搞到了12点半。发现坑数次。还专门用js写了个画圆的。当然我现在知道怎么画圆了。过两天去把oj7的统计图改一改。

思路是把所有圆的左右端点和所有圆的两两交的x提出来，unique一发（不然会tle显然的）。这样两个x之间的部分的圆弧是不会有交的。那么就可以视为一堆弓形和梯形了。然后扫描一下中位线的有效长度算梯形。每次新进入一个连通区域和出来的时候加相应弓形的面积。

需要注意按y排序的时候是按弓形与竖线的交点的y，而不是梯形的y，否则会有非常神奇的精度误差。另外算弓形面积必需要用反三角函数（至少我没研究出怎么不用）。如果用acos的话会暴精度暴得很惨。要用atan2来提高精度。

然后这样写连样例都会tle。仔细感受一下发现好像如果有一堆同心圆的话都会tle。然后考虑去掉包含的圆之后似乎就快了。然后就交了，然后居然跑过了。

然后来欣赏一下我的时间巨大，空间巨大，长度巨大。挤在一堆simpson里显得特别郁闷的代码。

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

using namespace std;

const double eps = 1e-8;
const double pi = acos(-1);

inline int sgn(double x) {
	return (x > eps) - (x < -eps);
}
inline double sqr(double x) {
	return x * x;
}

typedef struct geo_obj {
	double x, y;
	geo_obj() {}
	geo_obj(double xo, double yo) {
		x = xo, y = yo;
	}
	inline void read() {
		scanf("%lf%lf", &x, &y);
	}
	inline geo_obj vertical() {
		return geo_obj(-y, x);
	}
	inline double len() {
		return sqrt(sqr(x) + sqr(y));
	}
} point, vect;
inline geo_obj operator +(const geo_obj& a, const geo_obj& b) {
	return geo_obj(a. x + b. x, a. y + b. y);
}
inline geo_obj operator -(const geo_obj& a, const geo_obj& b) {
	return geo_obj(a. x - b. x, a. y - b. y);
}
inline geo_obj operator *(const geo_obj& a, const double& b) {
	return geo_obj(a. x * b, a. y * b);
}
inline double operator *(const geo_obj& a, const geo_obj& b) {
	return a. x * b. y - a. y * b. x;
}
inline double operator &(const geo_obj& a, const geo_obj& b) {
	return a. x * b. x + a. y * b. y;
}
inline double sqDist(const geo_obj& a, const geo_obj& b) {
	return sqr(a. x - b. x) + sqr(a. y - b. y);
}

struct circle {
	point o;
	double r, sr;
	inline void read() {
		o. read();
		scanf("%lf", &r);
		sr = sqr(r);
	}
};

inline double helenSqr(double a, double b, double c) {
	double p((a + b + c) / 2.0);
	return sqrt(p * (p - a) * (p - b) * (p - c));
}

struct ypr {
	double y, s, yv;
	int v;
};
inline bool operator <(const ypr& a, const ypr& b) {
	return (a. yv < b. yv) || (a. yv == b. yv && a. v > b. v);
}

const int maxn = 1009;
const int maxx = 1002009;

int n, totx;
double x[maxx], ans;
circle a[maxn];
ypr y[maxn << 1];

double arcSqr(circle c, point a, point b) {
	vect x(a - c. o), y(b - c. o);
	double sqt(fabs(x * y));
	double d(sqrt(sqDist(a, b)));
	double tht(atan2(d / 2, sqt / d) * 2);
	//double tht(acos((x & y) / c. sr));
	sqt /= 2;
	return tht * c. sr / 2 - sqt;
}

int main() {
#ifndef ONLINE_JUDGE
	freopen("in.txt", "r", stdin);
#endif

	scanf("%d", &n);
	totx = 0;
	for (int i = 0; i < n; ++ i) { 
		a[i]. read();
		bool thr(0);
		for (int j = 0; j < i && !thr; ++ j)
			if (sqDist(a[i]. o, a[j]. o) <= sqr(a[i]. r - a[j]. r)) {
				if (a[i]. r > a[j]. r)
					swap(a[i], a[j]);
				thr = 1;
			}
		if (thr) {
			-- i;
			-- n;
		}
	}
	for (int i = 0; i < n; ++ i) {
		x[totx ++] = a[i]. o. x - a[i]. r;
		x[totx ++] = a[i]. o. x + a[i]. r;
	}
	for (int i = 1; i < n; ++ i)
		for (int j = 0; j < i; ++ j) {
			double sd(sqDist(a[i]. o, a[j]. o));
			if (sd <= sqr(a[i]. r + a[j]. r) && sd >= sqr(a[i]. r - a[j]. r)) {
				double d(sqrt(sd));
				double s(helenSqr(a[i]. r, a[j]. r, d));
				double h(s * 2 / d);
				double l(sqrt(a[i]. sr - sqr(h)));
				if (sqr(a[j]. r) > sd + sqr(a[i]. r))
					l = -l;
				point p(a[i]. o + (a[j]. o - a[i]. o) * (l / d));
				vect v((a[j]. o - a[i]. o). vertical());
				v = v * (h / v. len());
				x[totx ++] = p. x + v. x;
				x[totx ++] = p. x - v. x;
			}
		}
	sort(x, x + totx);
	totx = unique(x, x + totx) - x;
	ans = 0;
	for (int i = 0; i < totx - 1; ++ i) {
		double xl(x[i]), xr(x[i + 1]);
		int cy(0);
		for (int j = 0; j < n; ++ j) 
			if (a[j]. o. x - a[j]. r < xr && a[j]. o. x + a[j]. r > xl) {
				double yl(sqrt(a[j]. sr - sqr(a[j]. o. x - xl)));
				double yv(sqrt(a[j]. sr - sqr(a[j]. o. x - (xl + xr) / 2.0)));
				double yr(sqrt(a[j]. sr - sqr(a[j]. o. x - xr)));
				double ar(arcSqr(a[j], point(xl, a[j]. o. y - yl), point(xr, a[j]. o. y - yr)));
				double ym((yl + yr) / 2);
				y[cy]. y = a[j]. o. y - ym;
				y[cy]. yv = a[j]. o. y - yv;
				y[cy]. v = 1;
				y[cy]. s = ar;
				++ cy;
				y[cy]. y = a[j]. o. y + ym;
				y[cy]. yv = a[j]. o. y + yv;
				y[cy]. v = -1;
				y[cy]. s = ar;
				++ cy;
			}
		double ml(0);
		sort(y, y + cy);
		for (int j = 0, cnt = 0; j < cy - 1; ++ j) {
			if ((y[j]. v == 1 && !cnt) || (y[j]. v == -1 && cnt == 1))
				ans += y[j]. s;
			cnt += y[j]. v;
			if (cnt)
				ml += y[j + 1]. y - y[j]. y;
		}
		ans += y[cy - 1]. s;
		ans += ml * (xr - xl);
	}
	printf("%.3lf\n", ans);
}