[BZOJ2178]圆的面积并
[BZOJ2178]圆的面积并
试题描述
给出 \(N\) 个圆，求其面积并
输入
先给一个数字 \(N\), \(N \le 1000\) 接下来是 \(N\) 行是圆的圆心，半径，其绝对值均为小于 \(1000\) 的整数
输出
面积并，保留三位小数
输入示例 & 输出示例
数据规模及约定
见“输入”
题解
裸的 Simpson 积分，可以用来解决一般的连续函数的积分。
Simpson 积分是这样一个式子
对于 \(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;
}
