【BZOJ 2178】圆的面积并 自适应辛普森公式

$\int_l^r {f(x) = \frac{{(r - l)(f(l) + f(r) + 4f(mid))}}{6}}$
f(i)是x=i这一条直线经过的被圆覆盖的长度，然后递归求解，如果发现误差小于1e-13就直接返回

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#define eps 1e-13
#define maxn 1021
using namespace std;
int n,can[maxn],cnt,vis[maxn];
struct P{
    double x,y;
    P(double a=0,double b=0):x(a),y(b){}
};
typedef P vec;
vec operator-(P a,P b){return vec(a.x-b.x,a.y-b.y);}
double operator*(vec a,vec b){return a.x*b.y-a.y*b.x;}
vec operator*(vec a,double t){return vec(a.x*t,a.y*t);}
double sqr(double x){return x*x;}
double dis(P a,P b){return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}
struct Cir{
    P o;double r;
    bool operator<(const Cir& b)const{
        return o.x-r < b.o.x-b.r;
    }
}c[maxn],a[maxn];
bool dcmp(double x,double y){return fabs(x-y)<eps;}
struct Line{
    double l,r;
    bool operator<(const Line& y)const{
        return dcmp(l,y.l) ? r<y.r : l<y.l;
    }
}l[maxn];
double f(double x){
    int tot=0;double ans=0;
    for(int i=1;i<=cnt;i++){
        if(a[i].o.x-a[i].r >=x || a[i].o.x+a[i].r<=x)continue;
        l[++tot].l=a[i].o.y-sqrt(sqr(a[i].r)-sqr(a[i].o.x-x));
        l[tot].r=a[i].o.y+sqrt(sqr(a[i].r)-sqr(a[i].o.x-x));
    }sort(l+1,l+1+tot);
    double L=l[1].l,R=l[1].r;
    for(int i=2;i<=tot;i++){
        if(l[i].l>R)ans+=R-L,L=l[i].l,R=l[i].r;
        else R=max(R,l[i].r);
    }
    ans+=R-L;
    return ans;
}
double calc(double len,double fl,double fr,double fmid){
    return len*(fl+fr+4*fmid)/6;
}
double simpson(double l,double r,double mid,double fl,double fr,double fmid,double s){
    double midl=(mid+l)/2,midr=(r+mid)/2;
    double f1=f(midl),f2=f(midr);
    double sl=calc(mid-l,fl,fmid,f1),sr=calc(r-mid,fmid,fr,f2);
    if(dcmp(s,sl+sr))return s;
    else return simpson(l,mid,midl,fl,fmid,f1,sl)+simpson(mid,r,midr,fmid,fr,f2,sr);
}

int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lf%lf%lf",&c[i].o.x,&c[i].o.y,&c[i].r);
    for(int i=1;i<=n;i++){
        bool ok=true;
        for(int j=1;j<=n;j++){
            if(i==j||vis[j])continue;
            if(dis(c[i].o,c[j].o)+c[i].r<c[j].r+eps){vis[i]=1,ok=false;break;}
        }
        if(ok)c[++cnt]=c[i];
    }n=cnt;double ans=0,fl,fr,fmid,mid;
    sort(c+1,c+1+n);
    for(int j,i=1;i<=n;i=j){
        double L=c[i].o.x-c[i].r,R=c[i].o.x+c[i].r;
        for(j=i+1;j<=n;j++){
            if(c[j].o.x-c[j].r>R+eps)break;
            R=max(R,c[j].o.x+c[j].r);
        }cnt=0;
        for(int k=i;k<j;k++)a[++cnt]=c[k];
        mid=(R+L)/2;
        fl=f(L),fr=f(R);
        fmid=f(mid);
        ans+=simpson(L,R,mid,fl,fr,fmid,calc(R-L,fl,fr,fmid));
    }
    printf("%.3lf",ans);
    return 0;
}