BZOJ2178: 圆的面积并(Simpson积分)

传送门

题意：
求圆的面积并。

题解：
Simpson积分。

先找出不同的块，对于一个块积分逼近。
用f(l,r)表示Simpson积分在区间[l,r]的值。
那么有：

S[l,r]={f(l,r)(fabs(f(l,mid)+f(mid,r)−f(l,r))<eps)f(l,mid)+f(mid,r)others

其中,f(l,r)=(len(l)+len(r)+4∗len(mid))∗(r−l)/6.0len(i)表示在横坐标i处直线与图像的并集的长度。

#include <bits/stdc++.h>using namespace std;const int Maxn = 1e3 + 50;const double eps = 1e-11;int read() {    char ch = getchar();    int i = 0, f = 1;    while (!isdigit(ch)) {        if (ch == '-') f = -1;        ch = getchar();    }    while (isdigit(ch)) {        i = (i << 1) + (i << 3) + ch - '0';        ch = getchar();    }    return i * f;}struct CIR {    double x, y, r;} cirtmp[Maxn], cir[Maxn];struct SEG {    double x, y;} seg[Maxn];int n, tot, del[Maxn], bg, ed;double ans;bool compr(const CIR &a, const CIR &b) { return a.r > b.r; }bool compx(const CIR &a, const CIR &b) { return a.x - a.r < b.x - b.r; }bool compseg(const SEG &a, const SEG &b) {    return (a.x == b.x) ? (a.y < b.y) : (a.x < b.x);}double dis(const CIR &a, const CIR &b) {    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));}void init() {    n = read();    for (int i = 1; i <= n; i++) {        cirtmp[i].x = read();        cirtmp[i].y = read();        cirtmp[i].r = read();    }    sort(cirtmp + 1, cirtmp + n + 1, compr);    for (int i = 1; i <= n; i++)        for (int j = i + 1; j <= n; j++)            if (dis(cirtmp[i], cirtmp[j]) <= cirtmp[i].r - cirtmp[j].r)                del[j] = 1;    for (int i = 1; i <= n; i++)        if (!del[i]) cir[++tot] = cirtmp[i];    sort(cir + 1, cir + tot + 1, compx);}bool calc(CIR a, double pos, SEG &SEGtmp) {    if (a.x + a.r < pos || a.x - a.r > pos) return false;    double t = a.r * a.r - (a.x - pos) * (a.x - pos);    if (t <= 0) return false;    double x1 = a.y - sqrt(t), x2 = a.y + sqrt(t);    SEGtmp.x = x1;    SEGtmp.y = x2;    return true;}double getlen(double pos) {    static int Segcnt, head, tail;    static SEG SEGtmp;    static double len;    Segcnt = 0;    len = 0;    for (int i = bg; i <= ed; i++)        if (calc(cir[i], pos, SEGtmp)) seg[++Segcnt] = SEGtmp;    if (!Segcnt) return 0;    sort(seg + 1, seg + Segcnt + 1, compseg);    double l, r;    for (int i = 1; i <= Segcnt; i = tail + 1) {        head = tail = i, l = seg[i].x, r = seg[i].y;        for (int j = i + 1; j <= Segcnt; j++) {            if (seg[j].x > r) break;            if (seg[j].y > r) r = seg[j].y;            tail++;        }        len += r - l;    }    return len;}double calc(double l, double r, double llen, double rlen, double midlen) {    return (llen + rlen + 4.0 * midlen) * (r - l) / 6.0;}double simpson(double l, double r, double S, double llen, double rlen,               double midlen) {    double mid = (l + r) / 2.0;    double lmid = (l + mid) / 2.0, rmid = (r + mid) / 2.0;    lmid = getlen(lmid), rmid = getlen(rmid);    double S1 = calc(l, mid, llen, midlen, lmid),           S2 = calc(mid, r, midlen, rlen, rmid);    if (fabs(S1 + S2 - S) < eps)        return S1 + S2;    else        return simpson(l, mid, S1, llen, midlen, lmid) +               simpson(mid, r, S2, midlen, rlen, rmid);}int main() {    init();    for (int i = 1; i <= tot; i = ed + 1) {        static double l, r;        bg = ed = i;        l = cir[i].x - cir[i].r;        r = cir[i].x + cir[i].r;        for (int j = i + 1; j <= tot; j++) {            if (cir[j].x - cir[j].r > r) {                break;            }            if (cir[j].x + cir[j].r > r) r = cir[j].x + cir[j].r;            ed = j;        }        double mid = (l + r) / 2;        double llen = getlen(l), rlen = getlen((r)), midlen = getlen(mid);        double nowS = calc(l, r, llen, rlen, midlen);        ans += simpson(l, r, nowS, llen, rlen, midlen);    }    printf("%0.3f", (double)ans);    return 0;}