# Area of Union of Circles

参考：https://codeforces.com/blog/entry/2111?#comment-488280

## 問題設定
円が $10^3$ 個くらい与えられるので和集合Dの面積を求めましょう。

## 考察
平面走査をしようにも途中で謎の図形が出てきて等積変形に困る。
そこで **Greenの定理** という飛び道具を用いることを考える。

> 多重連結領域 $D$ と向きづけられた境界 $C$ 、 $C^1$ 級関数 $P(x,y),Q(x,y)$ に対し、
> $\oint_C (Pdx+Qdy)=\int\int_D (\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}) dxdy$
> が成り立つ。

ここで $\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y}=1$ なる $P,Q$ (例えば $P=0,Q=x$ など) を取ると右辺は領域$D$ の面積に一致する。
よって $D$ の面積を求めることが左辺の線積分を実行することに言い換えられた。

Dの境界(反時計回りに向き付け)を所属している円ごとに分解すると、各々の円周のうち他の円の内部に入っていない円弧の和に一致することが分かる。(適当な図を書くと納得できると思います)
よってそのような円弧それぞれをパラメータ化して左辺の積分の和を計算すれば和集合の面積が計算できる。
計算量について、各円と他の円との交点を偏角ソートする部分がボトルネックになり$O(N^2\log N)$ になる。

(追記) $D$ の境界は $O(n)$ 個の円弧の和集合で表せるらしい(Power Diagram)ので、それを使うと $O(N\log N)$ になる。

## 実装

T UnionArea(vector<Circle> &cs) {
    int n = SZ(cs);
    double ret = 0;
    rep(i, 0, n) {
        vector<Point> pts;
        bool inside = 0;
        rep(j, 0, n) if (i != j) {
            if (eq(abs(cs[i].p - cs[j].p), 0.) and
                eq(abs(cs[i].r - cs[j].r), 0.)) {
                if (i > j) {
                    inside = 1;
                    break;
                }
            } else {
                double D = abs(cs[i].p - cs[j].p);
                if (cs[j].r + eps >= D + cs[i].r) {
                    inside = 1;
                    break;
                }
            }
            auto is = Intersection(cs[i], cs[j]);
            if (SZ(is) >= 2) {
                if (cross(is[0] - cs[i].p, cs[j].p - cs[i].p) < 0)
                    swap(is[0], is[1]);
                pts.push_back(is[0]);
                pts.push_back(is[1]);
            }
        }
        if (inside)
            continue;
        if (SZ(pts) == 0) {
            ret += pow(cs[i].r, 2) * M_PI;
            continue;
        }
        vector<int> ord(SZ(pts));
        iota(ALL(ord), 0);
        sort(ALL(ord), [&](int x, int y) {
            return cmp(pts[x] - cs[i].p, pts[y] - cs[i].p);
        });
        vector<int> rev(SZ(pts));
        rep(x, 0, SZ(pts)) rev[ord[x]] = x;
        vector<int> rui(SZ(pts) + 1);
        rep(x, 0, SZ(pts) / 2) {
            int st = rev[x * 2], en = rev[x * 2 + 1];
            if (st < en) {
                rui[st]++;
                rui[en]--;
            } else {
                rui[0]++;
                rui[en]--;
                rui[st]++;
            }
        }
        show(SZ(pts));
        rep(x, 0, SZ(pts)) rui[x + 1] += rui[x];
        rep(x, 0, SZ(pts)) {
            double theta = arg(pts[ord[x]] - cs[i].p);
        }
        rep(x, 0, SZ(pts)) if (rui[x] == 0) {
            double theta1 = arg(pts[ord[(x + 1) % SZ(pts)]] - cs[i].p);
            double theta2 = arg(pts[ord[x]] - cs[i].p);
            if (theta1 < theta2)
                theta1 += M_PI * 2;
            double add =
                cs[i].r * cs[i].p.X * sin(theta1) +
                pow(cs[i].r, 2) * (theta1 + sin(theta1) * cos(theta1)) / 2;
            add -= cs[i].r * cs[i].p.X * sin(theta2) +
                   pow(cs[i].r, 2) * (theta2 + sin(theta2) * cos(theta2)) / 2;
            ret += add;
        }
    }
    return ret;
}