bzoj 2178 圆的面积并——辛普森积分

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define db double
using namespace std;
const int N=1005; const db eps=1e-13;
db Fabs(db x){return x<0?-x:x;}
db Sqr(db x){return x*x;}
int nMx(int a,int b){return a>b?a:b;}
db dMx(db a,db b){return a>b?a:b;}
int n,xl[N],xr[N],st,ed; db ans;
struct Cir{
  int r,x,y;
  bool operator< (const Cir &b)const
  {return r>b.r;}
}t[N];
struct Node{
  db l,r;
  Node(db a=0,db b=0):l(a),r(b) {}
  bool operator< (const Node &b)const
  {return l<b.l;}
}a[N];
bool cmp(Cir u,Cir v){return u.x-u.r<v.x-v.r;}
db f(db p)
{
  int cnt=0; db ret=0;
  for(int i=st;i<=ed;i++)
    if(xr[i]>p&&xl[i]<p)
      {
    db k=sqrt(t[i].r-Sqr(p-t[i].x));
    a[++cnt]=Node(t[i].y-k,t[i].y+k);
      }
  sort(a+1,a+cnt+1);
  for(int i=1,j=1;i<=cnt;i=j)
    {
      db L=a[i].l,R=a[i].r;
      for(j=i+1;j<=cnt&&a[j].l<=R;j++)R=dMx(R,a[j].r);
      ret+=R-L;
    }
  return ret;
}
db cal(db fl,db fm,db fr,db d){return (fl+4*fm+fr)*d/6;}
db simp(db l,db mid,db r,db fl,db fm,db fr,db ret)
{
  db m1=(l+mid)*0.5,m2=(mid+r)*0.5,fm1=f(m1),fm2=f(m2);
  db vl=cal(fl,fm1,fm,mid-l),vr=cal(fm,fm2,fr,r-mid);
  if(Fabs(vl+vr-ret)<=eps)return vl+vr;
  return simp(l,m1,mid,fl,fm1,fm,vl)+simp(mid,m2,r,fm,fm2,fr,vr);
}
int main()
{
  scanf("%d",&n);
  for(int i=1;i<=n;i++)
    scanf("%d%d%d",&t[i].x,&t[i].y,&t[i].r);
  sort(t+1,t+n+1);
  for(int i=1;i<=n;i++)
    for(int j=1;j<i;j++)
      {
    db d=sqrt(Sqr(t[i].y-t[j].y)+Sqr(t[i].x-t[j].x));
    if(d+t[i].r<=t[j].r)
      {t[i--]=t[n--];break;}
      }
  sort(t+1,t+n+1,cmp);
  for(int i=1;i<=n;i++)
    xl[i]=t[i].x-t[i].r,xr[i]=t[i].x+t[i].r,t[i].r*=t[i].r;
  int L,R;db fl,fr,mid,fm;
  for(int i=1,j;i<=n;i=j)//j not j-1
    {
      L=xl[i];R=xr[i];
      for(j=i+1;j<=n&&xl[j]<=R;j++)R=nMx(R,xr[j]);
      st=i;ed=j-1;
      fl=f(L);fr=f(R);mid=(db)(L+R)*0.5;fm=f(mid);
      ans+=simp(L,mid,R,fl,fm,fr,cal(fl,fm,fr,R-L));
    }
  printf("%.3f
",ans);
  return 0;
}