BZOJ 2178 圆的面积并

自适应辛普森积分

别的不会，裸上辛普森

#include<cstdio>
#include<cmath>
#include<algorithm>
#define N 1005
using namespace std;
struct circle
{
    double r, x, y;
}c[N];
struct interval
{
    double l, r;
    friend bool operator < (interval a, interval b)
    {
        if(a.l!=b.l)return a.l<b.l;
        else return a.r>b.r;
    }
}inter[N];
const double INF = 1000000000.0, eps = 1e-13;
int n;
bool ban[N];
bool cover(circle a, circle b)
{
    double dis = sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
    return dis + b.r <= a.r;
}
bool cmp(circle a, circle b)
{
    if(a.x-a.r != b.x-b.r)return a.x-a.r < b.x-b.r;
    else return a.x+a.r > b.x+b.r;
}
double f(double x)
{
    int icnt=0;
    double d, h, ret=0, l = -INF;
    for(int i = 1; i <= n; i++)
    {
        d = abs(c[i].x-x);
        if(d>=c[i].r)continue;
        h = sqrt(c[i].r*c[i].r - d*d);
        inter[++icnt]=(interval){c[i].y-h, c[i].y+h};
    }
    if(!icnt)return 0;
    sort(inter+1,inter+1+icnt);
    for(int i = 1; i <= icnt; i++)
    {
        if(inter[i].r <= l)continue;
        l = max(l,inter[i].l);
        ret += inter[i].r-l;
        l = inter[i].r;
    }
    return ret;
}
double Simpson(double l, double r, double fl, double fm, double fr)
{
    double mid = (l+r)/2, lmid = (l+mid)/2, rmid = (mid+r)/2;
    double flm = f(lmid), frm = f(rmid);
    if(abs ((fl+4*flm+fm)/6*(mid-l) + (fm+4*frm+fr)/6*(r-mid) - (fl+4*fm+fr)/6*(r-l) ) < eps)
        return (fl+4*fm+fr)/6*(r-l);
    else return Simpson(l,mid,fl,flm,fm) + Simpson(mid,r,fm,frm,fr);
}
int main()
{
    scanf("%d",&n);
    for(int i = 1; i <= n; i++)
        scanf("%lf%lf%lf",&c[i].x,&c[i].y,&c[i].r);
    for(int i = 1; i <= n; i++)
        if(!ban[i])
            for(int j = i+1; j <= n; j++)
                if(!ban[j])
                    if(cover(c[i],c[j]))
                        ban[j]=1;
    int pos = 0;
    for(int i = 1; i <= n; i++)
        if(!ban[i])
            c[++pos]=c[i];
    n=pos;
    sort(c+1,c+1+n,cmp);
    double l=-INF, ans=0;
    for(int i = 1; i <= n; i++)
    {
        if(c[i].x + c[i].r <= l+eps)continue;
        l = max(l,c[i].x-c[i].r);
        ans += Simpson(l,c[i].x+c[i].r,f(l),f((l+c[i].x+c[i].r)/2),f(c[i].x+c[i].r));
        l = c[i].x + c[i].r;
    }
    printf("%.3lf\n",ans);
}