自适应辛普森积分

　　辛普森真是个强大的东西……很多东西都能积= =

　　这题的正解看上去很鬼畜，至少我这种不会计算几何的渣渣是写不出来……（对圆的交点求图包，ans=凸包的面积+一堆弓形的面积，另外还有中空的情况……那种凸包怎么求啊喂！）

  1 /**************************************************************
  2     Problem: 2178
  3     User: Tunix
  4     Language: C++
  5     Result: Accepted
  6     Time:8808 ms
  7     Memory:1372 kb
  8 ****************************************************************/
  9
 10 //BZOJ 2178
 11 #include<cmath>
 12 #include<vector>
 13 #include<cstdio>
 14 #include<cstring>
 15 #include<cstdlib>
 16 #include<iostream>
 17 #include<algorithm>
 18 #define rep(i,n) for(int i=0;i<n;++i)
 19 #define F(i,j,n) for(int i=j;i<=n;++i)
 20 #define D(i,j,n) for(int i=j;i>=n;--i)
 21 #define pb push_back
 22 using namespace std;
 23 inline int getint(){
 24     int v=0,sign=1; char ch=getchar();
 25     while(ch<‘0‘||ch>‘9‘){ if (ch==‘-‘) sign=-1; ch=getchar();}
 26     while(ch>=‘0‘&&ch<=‘9‘){ v=v*10+ch-‘0‘; ch=getchar();}
 27     return v*sign;
 28 }
 29 const int N=1010,INF=~0u>>2;
 30 const double eps=1e-13;
 31 typedef long long LL;
 32 typedef double db;
 33 /******************tamplate*********************/
 34 int n,top,st,ed;
 35 db xl[N],xr[N],ans;
 36 bool del[N];
 37 struct data{db x,y,r;}t[N],sk[N];
 38 struct line{db l,r;}p[N];
 39 db dis(data a,data b){
 40     return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
 41 }
 42 bool cmp1(data a,data b){return a.r<b.r;}
 43 bool cmp2(data a,data b){return a.x-a.r<b.x-b.r;}
 44 bool cmp3(line a,line b){return a.l<b.l;}
 45 void init(){
 46     scanf("%d",&n);
 47     F(i,1,n) scanf("%lf%lf%lf",&t[i].x,&t[i].y,&t[i].r);
 48     sort(t+1,t+n+1,cmp1);
 49     F(i,1,n) F(j,i+1,n)
 50         if(dis(t[i],t[j])<=t[j].r-t[i].r){
 51             del[i]=1; break;
 52         }
 53     F(i,1,n)if(!del[i]) sk[++top]=t[i];n=top;
 54     sort(sk+1,sk+n+1,cmp2);
 55 }
 56 db getf(db x){
 57     int j,sz=0; db r,len=0,dis;
 58     F(i,st,ed){
 59         if (x<=xl[i]||x>=xr[i])continue;
 60         dis=sqrt(sk[i].r-(x-sk[i].x)*(x-sk[i].x));
 61         p[++sz].l=sk[i].y-dis; p[sz].r=sk[i].y+dis;
 62     }
 63     sort(p+1,p+sz+1,cmp3);
 64     F(i,1,sz){
 65         r=p[i].r;
 66         for(j=i+1;j<=sz;j++){
 67             if(p[j].l>r)break;
 68             if (r<p[j].r) r=p[j].r;
 69         }
 70         len+=r-p[i].l; i=j-1;
 71     }
 72     return len;
 73 }
 74 db cal(db l,db fl,db fmid,db fr){
 75     return (fl+fmid*4+fr)*l/6;
 76 }
 77 db simpson(db l,db mid,db r,db fl,db fmid,db fr,db s){
 78     db m1=(l+mid)/2,m2=(mid+r)/2;
 79     db f1=getf(m1),f2=getf(m2);
 80     db g1=cal(mid-l,fl,f1,fmid),g2=cal(r-mid,fmid,f2,fr);
 81     if (fabs(g1+g2-s)<eps)return g1+g2;
 82     return simpson(l,m1,mid,fl,f1,fmid,g1)+
 83            simpson(mid,m2,r,fmid,f2,fr,g2);
 84 }
 85 void work(){
 86     int i,j;
 87     db l,r,mid,fl,fr,fmid;
 88     F(i,1,n){
 89         xl[i]=sk[i].x-sk[i].r;
 90         xr[i]=sk[i].x+sk[i].r;
 91         sk[i].r*=sk[i].r;
 92     }
 93     F(i,1,n){
 94         l=xl[i]; r=xr[i];
 95         for(j=i+1;j<=n;j++){
 96             if (xl[j]>r) break;
 97             if (xr[j]>r) r=xr[j];
 98         }
 99         st=i; ed=j-1;i=j-1;
100         mid=(l+r)/2;
101         fl=getf(l); fr=getf(r); fmid=getf(mid);
102         ans+=simpson(l,mid,r,fl,fmid,fr,cal(r-l,fl,fmid,fr));
103     }
104 }
105 int main(){
106 #ifndef ONLINE_JUDGE
107     freopen("2178.in","r",stdin);
108     freopen("2178.out","w",stdout);
109 #endif
110     init();
111     work();
112     printf("%.3lf",ans);
113     return 0;
114 }