BZOJ 2178: 圆的面积并 [辛普森积分 区间并]

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=1005,INF=1e9;
const double eps=1e-13;
inline int sgn(double x){
    if(abs(x)<eps) return 0;
    else return x<0?-1:1;
}
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1; c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘; c=getchar();}
    return x*f;
}
int n;
struct Interval{
    double l,r;
    bool operator <(const Interval &a) const{
        return l==a.l?r<a.r:l<a.l;
    }
    Interval(double a=0,double b=0):l(a),r(b){}
}a[N],b[N];

struct Circle{
    int x,y,r;
    Circle(){}
    Circle(int x,int y,int r):x(x),y(y),r(r){}
    Interval f(double p){
        double dx=x-p,t=sqrt(r*r-dx*dx);
        return Interval(y-t,y+t);
    }
}C[N];
double Dis(Circle a,Circle b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool del[N];
void iniCircle(){
    for(int i=1;i<=n;i++) if(!del[i])
        for(int j=i+1;j<=n;j++) if(!del[j])
            if(Dis(C[i],C[j])<=abs(C[i].r-C[j].r)){
                if(C[i].r>C[j].r) del[j]=1;
                else del[i]=1;
            }
    int p=0;
    for(int i=1;i<=n;i++) if(!del[i]) C[++p]=C[i];
    n=p;
}
inline double F(double x){
    int m=0;
    for(int i=1;i<=n;i++)
        if(C[i].x-C[i].r<=x&&x<=C[i].x+C[i].r) a[++m]=C[i].f(x);
    sort(a+1,a+1+m);
    double last=-INF,re=0;
    for(int i=1;i<=m;i++){
        if(a[i].l>last) re+=a[i].r-a[i].l,last=a[i].r;
        else if(a[i].r>last) re+=a[i].r-last,last=a[i].r;
    }
    return re;
}

inline double cal(double l,double r){
    return (F(l)+F(r)+4*F((l+r)/2))*(r-l)/6;
}
inline double cal(double l,double r,double fl,double fr,double fm){
    return (fl+fr+4*fm)*(r-l)/6;
}
double Simpson(double l,double r,double now,double fl,double fr,double fm){
    double mid=(l+r)/2,flm=F((l+mid)/2),frm=F((mid+r)/2);
    double p=cal(l,mid,fl,fm,flm),q=cal(mid,r,fm,fr,frm);
    if(sgn(now-p-q)==0) return now;
    else return Simpson(l,mid,p,fl,fm,flm)+Simpson(mid,r,q,fm,fr,frm);
}

double lb=INF,rb=-INF;
int main(int argc, const char * argv[]) {
    n=read();
    for(int i=1;i<=n;i++){//printf("i %d\n",i);
        C[i].x=read();C[i].y=read();C[i].r=read();
        lb=min(lb,(double)C[i].x-C[i].r);
        rb=max(rb,(double)C[i].x+C[i].r);
    }
    iniCircle();
    //printf("%d\n",n);
    double fl=F(lb),fr=F(rb),fm=F((lb+rb)/2);
    printf("%.3lf",Simpson(lb,rb,cal(lb,rb,fl,fr,fm),fl,fr,fm));
    return 0;
}