BZOJ2178 圆的面积并

			    <div class="post_brief"><p>
计算几何题。别人都用simpson搞。然后我决定不水一发，用hwd神犇在wc上讲的关键点法做。</p>

然后就搞到了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", &amp;n);
totx = 0;
for (int i = 0; i &lt; n; ++ i) { 
	a[i]. read();
	bool thr(0);
	for (int j = 0; j &lt; i &amp;&amp; !thr; ++ j)
		if (sqDist(a[i]. o, a[j]. o) &lt;= sqr(a[i]. r - a[j]. r)) {
			if (a[i]. r &gt; a[j]. r)
				swap(a[i], a[j]);
			thr = 1;
		}
	if (thr) {
		-- i;
		-- n;
	}
}
for (int i = 0; i &lt; 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 &lt; n; ++ i)
	for (int j = 0; j &lt; i; ++ j) {
		double sd(sqDist(a[i]. o, a[j]. o));
		if (sd &lt;= sqr(a[i]. r + a[j]. r) &amp;&amp; sd &gt;= 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) &gt; 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 &lt; totx - 1; ++ i) {
	double xl(x[i]), xr(x[i + 1]);
	int cy(0);
	for (int j = 0; j &lt; n; ++ j) 
		if (a[j]. o. x - a[j]. r &lt; xr &amp;&amp; a[j]. o. x + a[j]. r &gt; 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 &lt; cy - 1; ++ j) {
		if ((y[j]. v == 1 &amp;&amp; !cnt) || (y[j]. v == -1 &amp;&amp; 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);

}