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

∫rlf(x)=(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 1021using 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;}