bzoj2178 圆的面积并（自适应Simpson积分）

bzoj2178 圆的面积并

原题地址：http://www.lydsy.com/JudgeOnline/problem.php?id=2178

题意：
给出N个圆，求其面积并

数据范围
N< = 1000 ，圆心，半径，其绝对值均为小于1000的整数

题解：
当然simpson积分是很容易卡的，但是此题居然可过。
不过精度要开到1e-13。
一定要把被完全包含的圆去掉，这样可以快很多。

本来应该跳过空白分段去做，否则只要首位中都为空就直接返回0了，但此题坐标范围较小，没有卡。

代码：

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#define Poi Vec
using namespace std;
const double EPS=1e-13;
const int N=1005;
struct Vec
{
    double x,y;
    Vec(){}
    Vec(double x,double y):x(x),y(y){}
}L[N];
struct Cir
{
    Poi p;
    double r;
}C[N],O[N];
bool cmp_sig(const Vec &A,const Vec &B)
{return (A.x!=B.x)?A.x<B.x:A.y<B.y;}
bool cmp_cir(const Cir &A,const Cir &B) 
{return (A.p.x-A.r!=B.p.x-B.r)?A.p.x-A.r<B.p.x-B.r:A.p.x+A.r>B.p.x+B.r;}
double dis(const Poi &A,const Poi &B) {return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));}
bool del[N];
int n,top;
int sign(double x)
{
    if(fabs(x)<=EPS) return 0;
    else return x<0?-1:1;
}
double F(double x)
{
    top=0; double ret=0;
    for(int i=1;i<=n;i++)
    {
        double d=fabs(C[i].p.x-x);
        if(sign(d-C[i].r)>=0) continue;
        double s=sqrt(C[i].r*C[i].r-d*d);
        L[++top]=Vec(C[i].p.y-s,C[i].p.y+s);
    }
    if(!top) return 0;
    sort(L+1,L+top+1,cmp_sig);
    double l=L[1].x; double r=L[1].y;
    for(int i=2;i<=top;i++)
    {
        if(L[i].x<r) r=max(r,L[i].y);
        else {ret+=r-l; l=L[i].x;r=L[i].y;}
    }
    ret+=r-l;
    return ret;
}

double simpson(double a,double c,double b,double fa,double fc,double fb)
{return (fa+4*fc+fb)*(b-a)/6;}
double adjust(double a,double c,double b,double fa,double fc,double fb)
{
    double e=(double)(a+c)/2; double fe=F(e); 
    double f=(double)(c+b)/2; double ff=F(f);
    double lf=simpson(a,e,c,fa,fe,fc); double rg=simpson(c,f,b,fc,ff,fb);
    double alt=simpson(a,c,b,fa,fc,fb);
    if(sign(lf+rg-alt)==0) return alt;
    else return adjust(a,e,c,fa,fe,fc)+adjust(c,f,b,fc,ff,fb);
}
int main()
{
    scanf("%d",&n) ;double lf,rg; lf=2005.0; rg=-2005.0;
    for(int i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&O[i].p.x,&O[i].p.y,&O[i].r);
        lf=min(lf,O[i].p.x-O[i].r); rg=max(rg,O[i].p.x+O[i].r);
    }
    sort(O+1,O+n+1,cmp_cir);
    for(int i=1;i<=n;i++) 
    {   
        if(del[i]) continue;
        for(int j=i+1;j<=n;j++)
        {
            if(del[j]) continue;
            if(sign(dis(O[i].p,O[j].p)+O[j].r-O[i].r)<=0) del[j]=1;
        }
    }
    top=0; for(int i=1;i<=n;i++) if(!del[i]) C[++top]=O[i]; n=top;  
    double mid=(lf+rg)/2.0; double fa=F(lf),fb=F(rg),fc=F(mid);
    printf("%0.3lf\n",adjust(lf,mid,rg,fa,fc,fb));
    return 0;
}