BZOJ 2178: 圆的面积并 (辛普森积分)

code

#include <set>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int MAXN = 1005;
const double Pi = acos(-1.0);
const double eps = 1e-12;
const double INF = 1e12;
inline int dcmp(double x) {
	if(x <= eps && x >= -eps) return 0;
	return x > 0 ? 1 : -1;
}
inline double sqr(double x) { return x*x; }
struct Vector {
	double x, y;
	inline Vector(double _x=0, double _y=0): x(_x), y(_y){}
	inline Vector operator *(const double &k)const { return Vector(x*k, y*k) ;}
	inline Vector operator +(const Vector &t)const { return Vector(x+t.x, y+t.y); }
	inline Vector operator -(const Vector &t)const { return Vector(x-t.x, y-t.y); }
};
struct Point {
    double x, y;
    inline Point(double _x=0, double _y=0): x(_x), y(_y){}
    inline Vector operator -(const Point &t)const {  return Vector(x - t.x, y - t.y); }
    inline Point operator +(const Vector &t)const {  return Point(x + t.x, y + t.y); }
    inline double dist(const Point &t) { return sqrt(sqr(x-t.x) + sqr(y-t.y)); }
    inline double dist(const double &a, const double &b) { return sqrt(sqr(x-a) + sqr(y-b)); }
};
struct Line {
	Point p;
	Vector v;
	double ang;
	inline Line(Point _p=Point(0, 0), Vector _v=Vector(0, 0)): p(_p), v(_v), ang(atan2(v.y, v.x)){}
	inline bool operator <(const Line &t)const { return ang < t.ang; }
};
struct Circle {
	Point o; double r;
	inline Circle(Point _o=Point(0, 0), double _r=0): o(_o), r(_r){}
};
inline double Cross(const Vector &a, const Vector &b) { return a.x*b.y - a.y*b.x; }
inline bool Turn_Left(const Point &a, const Point &b, const Point &c) { return dcmp(Cross(b-a, c-a)) > 0; }
inline bool On_Left(const Line &a, const Point &b) { return dcmp(Cross(a.v, b-a.p)) >= 0; }
inline Point GLI(const Point &P, const Vector &v, const Point &Q, const Vector &w) { //Intersection 
	Vector u = Q - P;
	double k = Cross(u, w) / Cross(v, w);
	return P + v*k;
}
inline double Angle_C(const double &a, const double &b, const double &c) { //Triangle's 3 edges
	return acos((sqr(a)+sqr(b)-sqr(c))/(2*a*b));
}
inline double Tri_S(const Point &a, const Point &b, const Point &c) { return fabs(Cross(b-a, c-a))/2; }
inline Point rotate(const double &x, const double &y, const double °ree) {
	return Point(x*cos(degree)-y*sin(degree), x*sin(degree)+y*cos(degree));
}
int n; Circle cir[MAXN], tmp[MAXN];
inline bool cmp(const Circle &a, const Circle &b) {
	return a.o.x-a.r < b.o.x-b.r;
	//stO     Orz  -.-  >.<  xp   XD   qwq  QAQ   ToT  QwQ  QuQ   :-)   )-:   o-:   qx  TnT
}
inline bool cmp2(const Circle &a, const Circle &b) {
	return a.r > b.r;
}
inline bool chk(const Circle &a, const Circle &b) {
	return sqr(a.o.x-b.o.x) + sqr(a.o.y-b.o.y) <= sqr(a.r-b.r);
}
inline void Pre_work() {
	sort(cir + 1, cir + n + 1, cmp2);
	int cnt = 0;
	for(int i = 1; i <= n; ++i) if(cir[i].r) {
		bool flg = 1;
		for(int j = 1; j <= cnt; ++j)
			if(chk(cir[i], tmp[j])) { flg = 0; break; }
		if(flg) tmp[++cnt] = cir[i];
	}
	n = cnt;
	for(int i = 1; i <= n; ++i) cir[i] = tmp[i];
	sort(cir + 1, cir + n + 1, cmp);
}
int st, ed, tot;
pair<double, double>intervals[MAXN<<1];

inline void getCircleIntersection(const Circle &O, const double &x) {
	double len = sqrt(sqr(O.r) - sqr(x-O.o.x));
	intervals[++tot] = make_pair(O.o.y-len, O.o.y+len);
}

inline double f(double x) {
	tot = 0;
	for(int i = st; i <= ed; ++i)
		if(x < cir[i].o.x+cir[i].r && x > cir[i].o.x-cir[i].r)
			getCircleIntersection(cir[i], x);
	sort(intervals + 1, intervals + tot + 1);
	double L = -INF, R = -INF, res = 0;
	for(int i = 1; i <= tot; ++i)
		if(intervals[i].first > R)
			res += R - L, L = intervals[i].first, R = intervals[i].second;
		else R = max(R, intervals[i].second);
	return res + (R - L);
}
inline double Simpson(double L, double M, double R, double fL, double fM, double fR, int dep) {
    double M1 = (L + M) / 2, M2 = (M + R) / 2;
    double fM1 = f(M1), fM2 = f(M2);
    double g1 = (M-L) * (fL + 4*fM1 + fM) / 6, 
        g2 = (R-M) * (fM + 4*fM2 + fR) / 6,
        g = (R-L) * (fL + 4*fM + fR) / 6;
    if(dep > 11 && fabs(g-g1-g2) < 1e-10) return g; 
    else return Simpson(L, M1, M, fL, fM1, fM, dep+1) + Simpson(M, M2, R, fM, fM2, fR, dep+1);
}
inline double solve() {
	double ans = 0;
	for(int i = 1, j; i <= n; ++i) {
		double L = cir[i].o.x - cir[i].r, R = cir[i].o.x + cir[i].r;
		for(j = i+1; j <= n; ++j)
			if(cir[j].o.x - cir[j].r > R) break;
			else R = max(R, cir[j].o.x + cir[j].r);
		double M = (L + R) / 2;
		st = i, ed = i = j-1;
		double fL = f(L), fM = f(M), fR = f(R);
		ans += Simpson(L, M, R, fL, fM, fR, 0);
	}
	return ans;
}
int main () {
	scanf("%d", &n);
	for(int i = 1; i <= n; i++)
		scanf("%lf%lf%lf", &cir[i].o.x, &cir[i].o.y, &cir[i].r);
	Pre_work();
	printf("%.3f\n", solve());
}