uoj419 - 圆形 题解

考虑非辛普森积分的求圆的面积并的方法。毕竟这题就算不要求动态加入，\(n\leq 2000,\epsilon=10^{-9}\) 辛普森也别想跑了。

对于一个连通区域，如果在其中任作一条封闭曲线，其内部都属于这个区域，则称其为单连通区域。对应地，如果里面有「洞」，则称为复连通区域。

对于单连通区域求面积，有一个经典的方法：格林公式。具体公式涉及到我看不懂的微积分符号，这里仅用通俗的语言描述。

单连通区域显然有且仅有一条不自交的封闭曲线作为轮廓，设其为 \(L\)。任选一点 \(O\)，令动点 \(A\) 在 \(L\) 上逆时针移动一周，线段 \(OA\) 扫过的有向面积即为原图形的面积。有向面积的计算法则是什么？考虑 \(A\) 在某个具体位置 \(A_0\) 时，令 \(f(A_0)\) 的绝对值为 \(OA_0\) 的距离，符号取决于在移动 \(\epsilon\) 距离到达 \(A_1\) 后，\(OA_1\) 相比 \(OA_0\) 是否在逆时针方向，最后将 \(f(A_0)\) 在 \(L\) 上积起来就是答案。

证明不会，因为我根本不会微积分，笑死。可以发现一般多边形求面积方法就是格林公式的一个特例。

可以看到，格林公式的好处在于将一个单连通区域的面积转化为通过其轮廓就能计算出的信息，且易分解：必要时可以选择将轮廓断成几段，计算 \(A\) 经过每段扫过的有向面积相加。

回到原题，如果圆们的并是由若干单连通区域组成的话，那做法就显而易见了：考虑一个圆的弧的哪些部分暴露在外面（也就是作为轮廓）。只要枚举其它圆，看能覆盖哪段弧（这个在 P4518 里有，不再赘述），最后区间并一下，最后得到的区间相当于是要删去贡献的。为了方便，不妨先计算所有圆周都作为轮廓时的贡献（此时显然就是所有圆的面积和），然后删去该删的圆弧。

如何计算一段圆弧的贡献（扫过的有向面积）？考虑两端 \(X,Y\)（逆时针顺序），容易发现由一个面积一定是正的弓形 \(X\overset{\LARGE\frown}{YX}\)，加上 \(\triangle OXY\) 的有向面积组成，是 \(\mathrm O(1)\) 的。

如果连通区域内有洞怎么办呢？容易发现上述方法仍然是正确的：对一个洞，其轮廓的每段弧按其所在圆逆时针方向算贡献的话，它对于洞本身其实是顺时针方向的，恰好把这个洞的面积减掉了。

动态加入圆怎么做？很简单，只要用 set 随便维护每个圆的删除区间即可。一个需要注意的细节：若有两个圆重合，那么关于圆弧是否暴露在外面的判定会出 bug，需要手动 ban 掉重复圆，而其他情况都没问题，因为两个圆如果有大于 \(0\) 的重合长度，那这两个圆一定重合。

code

constexpr int N = 2010;
const db pi = acos(-1);

int sq(int x) { return x * x; }
db sq(db x) { return x * x; }

int n;
int xx[N], yy[N], ra[N];
bool ban[N];

set<tuple<db, db>> st[N];

void mian() {
  read();
  n = read();
  REP(i, 1, n) xx[i] = read(), yy[i] = read(), ra[i] = read();
  db ans = 0;
  auto deal = [&](int i, int coef, db l, db r) {
    ans -= coef * (r - l) / 2 * ra[i] * ra[i];
    ans += coef * ra[i] * ra[i] * sin(r - l) / 2;
    db x1 = xx[i] + ra[i] * cos(l), y1 = yy[i] + ra[i] * sin(l);
    db x2 = xx[i] + ra[i] * cos(r), y2 = yy[i] + ra[i] * sin(r);
    db cross = x1 * y2 - x2 * y1;
    ans -= coef * cross / 2;
  };
  auto ins = [&](int i, db l, db r) {
    auto &s = st[i];
    auto fd = s.lower_bound(mt(l, -1));
    while(fd != s.end() && X(*fd) <= r) {
      deal(i, -1, X(*fd), Y(*fd));
      chkmx(r, Y(*fd));
      s.erase(fd++);
    }
    if(fd != s.begin()) {
      --fd;
      if(Y(*fd) >= r) return;
      if(Y(*fd) >= l) {
        deal(i, -1, X(*fd), Y(*fd));
        l = X(*fd);
        s.erase(fd);
      }
    }
    deal(i, 1, l, r);
    s.insert(mt(l, r));
  };
  auto add = [&](int i, int j) {
    db dist = sqrt(sq(xx[i] - xx[j]) + sq(yy[i] - yy[j]));
    db sh = abs(dist - ra[i]), lo = dist + ra[i];
    if(lo <= ra[j]) return ins(i, 0, 2 * pi);
    if(sh >= ra[j]) return;
    db alpha = acos((sq(ra[i]) + sq(dist) - sq(ra[j])) / 2 / ra[i] / dist);
    db ang = atan2(yy[j] - yy[i], xx[j] - xx[i]);
    db l = ang - alpha, r = ang + alpha;
    while(l < 0) l += 2 * pi, r += 2 * pi;
    if(r >= 2 * pi) ins(i, l, 2 * pi), ins(i, 0, r - 2 * pi);
    else ins(i, l, r);
  };
  REP(i, 1, n) {
    REP(j, 1, i - 1) if(tie(xx[i], yy[i], ra[i]) == tie(xx[j], yy[j], ra[j])) ban[i] = true;
    if(!ban[i]) {
      ans += pi * ra[i] * ra[i];
      REP(j, 1, i - 1) if(!ban[j]) {
        add(i, j), add(j, i);
      }
    }
    printf("%.10Lf\n", ans);
  }
}