[BZOJ2178]圆的面积并

[BZOJ2178]圆的面积并

试题描述

给出 \(N\) 个圆，求其面积并

输入

先给一个数字 \(N\), \(N \le 1000\) 接下来是 \(N\) 行是圆的圆心，半径，其绝对值均为小于 \(1000\) 的整数

输出

面积并，保留三位小数

输入示例 & 输出示例

略

数据规模及约定

见“输入”

题解

裸的 Simpson 积分，可以用来解决一般的连续函数的积分。

Simpson 积分是这样一个式子

\[\int_l^r f(x) \mathrm{d} x \approx \frac{1}{6}(r-l)(f(l)+f(r)+4f(\frac{l+r}{2})) \]

对于 \(2\) 次及以下的函数上面的式子求出来的是精确值。

无聊的话可以证一下，令 \(f(x) = ax^2 + bx + c\)

\begin{aligned}
\int_l^r f(x) \mathrm{d} x
& = \int_l^r (ax^2 + bx + c) \mathrm{d}x \\
& = \left( \frac{1}{3} ar^3 + \frac{1}{2} br^2 + cr \right) - \left( \frac{1}{3} al^3 + \frac{1}{2} bl^2 + cl \right) \\
& = \frac{1}{3} a (r-l)(r^2 + lr + l^2) + \frac{1}{2} b (r-l)(r+l) + c(r-l) \\
& = \frac{1}{6} (r-l) ( 2ar^2 + 2alr + 2al^2 + 3br + 3bl + 6c ) \\
& = \frac{1}{6} (r-l) \left[ (ar^2 + br + c) + (al^2 + bl + c) + 4a \left( \frac{r}{2} \right)^2 + 4a\left( \frac{l}{2} \right)^2 + 4b\left( \frac{r}{2} \right) + 4b\left( \frac{l}{2} \right) + 4c \right] \\
& = \frac{1}{6} (r-l) \left[ f(r) + f(l) + 4f \left( \frac{l+r}{2} \right) \right]
\end{aligned}

对于这道题，想象一条竖直扫描线从左到右扫过整个图形，那么答案其实就是对“直线在图形内部的长度”进行积分，套用上面的公式就可以得出大致的答案了；为了提高精度，我们可以进行分治，对于横坐标在 \([l, r]\) 内的部分计算一下 \([l, r]\) 的答案，看和 \([l, m]\) 与 \([m, r]\)（\(m = \frac{l + r}{2}\)）的答案加起来是否相差很小，如果不够小则继续分治直接求 \([l, m]\) 与 \([m, r]\) 的和。

如何求出“在图形内部的长度”呢？注意到每个圆能转化成这条直线上的一个区间，最后就相当于求几个区间的并。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <cmath>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)
 
int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}
 
#define maxn 1010
 
const double eps = 1e-13;
struct Vec {
	double x, y;
	Vec() {}
	Vec(double _, double __): x(_), y(__) {}
	Vec operator - (const Vec& t) const { return Vec(x - t.x, y - t.y); }
	double len() { return sqrt(x * x + y * y); }
};
struct Circle {
	Vec p; double r;
	Circle() {}
	Circle(double _1, double _2, double _3): p(Vec(_1, _2)), r(_3) {}
} inc[maxn], cs[maxn];
int n;
bool del[maxn];
 
struct Line {
	double l, r;
	Line() {}
	Line(double _, double __): l(_), r(__) {}
	bool operator < (const Line& t) const { return l != t.l ? l < t.l : r > t.r; }
} ls[maxn];
int cl;
 
double calcL(double x) {
	cl = 0;
	rep(i, 1, n) {
		double d = abs(x - cs[i].p.x);
		if(d >= cs[i].r) continue;
		double h = sqrt(cs[i].r * cs[i].r - d * d);
		ls[++cl] = Line(cs[i].p.y - h, cs[i].p.y + h);
	}
	if(!cl) return 0.0;
	sort(ls + 1, ls + cl + 1);
	double nl = ls[1].l, nr = ls[1].r, sum = 0;
	rep(i, 2, cl) if(ls[i].l > nr) {
		sum += nr - nl;
		nl = ls[i].l; nr = ls[i].r;
	}
	else nr = max(nr, ls[i].r);
	return sum + nr - nl;
}
double solve(double l, double r, double lLen, double mLen, double rLen) {
	double mid = (l + r) * 0.5, lmLen = calcL((l + mid) * 0.5), rmLen = calcL((mid + r) * 0.5), S, lS, rS;
	S = (r - l) * (lLen + rLen + 4.0 * mLen) / 6.0;
	lS = (mid - l) * (lLen + mLen + 4.0 * lmLen) / 6.0;
	rS = (r - mid) * (mLen + rLen + 4.0 * rmLen) / 6.0;
	if(abs(lS + rS - S) < eps) return S;
	return solve(l, mid, lLen, lmLen, mLen) + solve(mid, r, mLen, rmLen, rLen);
}
 
int main() {
	double L = 10000, R = -10000;
	n = read();
	rep(i, 1, n) {
		double x = read(), y = read(), r = read();
		inc[i] = Circle(x, y, r);
		L = min(L, x - r); R = max(R, x + r);
	}
	rep(i, 1, n)
		rep(j, 1, n) if(inc[i].r > inc[j].r && (inc[i].p - inc[j].p).len() + inc[j].r <= inc[i].r) del[j] = 1;
	int nn = 0;
	rep(i, 1, n) if(!del[i]) cs[++nn] = inc[i];
	n = nn;
	 
	printf("%.3lf\n", solve(L, R, 0, calcL((L + R) * 0.5), 0));
	 
	return 0;
}