SPOJ - CIRU,bzoj 2178 圆的面积并

SPOJ-CIRU,bzoj2178 圆的面积并

题意

给出n个圆,求出它们的面积并.

做法

首先要引入辛普森积分:

用于求解由平滑曲线构成的面积,f(x)根据具体题意而定.
在本题中a,b表示求的组合图形的左右两端的横坐标,f(i)是x=i这条直线经过的图形的总长度.
比如说要用这个方法求一个圆心坐标为(1,2),半径为1的圆:
那么a=x-r=0,b=x+r=2,f(a)=0,f(b)=0,f((a+b)/2)=2*r=2.
所以原式=(b-a)/6*(f(a)+4 *f((a+b)/2)+f(b))=8/3,约为2.67
可以发现,这样直接算的精度很低,且(b-a)越大,就越不准确,因此我们在实际操作时,可以先算出a~b这段的面积并,再用这种方法分别算出a~(a+b)/2和(a+b)/2~b两段面积的和,比较两次计算的答案,若差值小于eps(自己定),则认为精度够了,返回,否则分别计算a~(a+b)/2和(a+b)/2~b两段面积的和
首先去掉被包含的圆,然后根据横坐标(圆不一定要相交,但横坐标一定要有交集)将组合图形分成很多部分,每个部分再用上述方法求出面积即可.

代码

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<vector>
#define db double
#define P pair<db,db>
#define mp make_pair
#define fi first
#define se second
#define N 1010
#define eps 1e-12
using namespace std;

int n,nn,s,t;
db xl[N],xr[N],ans;
bool xc[N];
struct Rd
{
    db x,y,r;
    bool operator < (const Rd u) const
    {
        return x-r<u.x-u.r;
    }
}rd[N],tmp[N];
vector<P>vec;

inline db getf(db x)
{
    int i,j;
    db tmp,res=0,l,r;
    vec.clear();
    for(i=s;i<=t;i++)
    {
        if(fabs(rd[i].x-x)>=rd[i].r) continue;
        tmp=sqrt(rd[i].r*rd[i].r-(rd[i].x-x)*(rd[i].x-x));
        vec.push_back(mp(rd[i].y-tmp,rd[i].y+tmp));
    }
    sort(vec.begin(),vec.end());
    for(i=0;i<vec.size();)
    {
        l=vec[i].fi,r=vec[i].se;
        for(j=i+1;j<vec.size();j++)
        {
            if(vec[j].fi>r) break;
            r=max(r,vec[j].se);
        }
        i=j;
        res+=r-l;
    }
    return res;
}

inline db len(Rd u,Rd v){return sqrt((u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y));}
inline db sim(db l,db r,db a,db b,db c){return (r-l)/6.0*(a+4.0*b+c);}

inline db calc(db l,db r,db mid,db fl,db fr,db fmid,db la)
{
    db a1,a2,m1,m2,fm1,fm2;
    m1=(l+mid)/2.0,m2=(mid+r)/2.0;
    fm1=getf(m1),fm2=getf(m2);
    a1=sim(l,mid,fl,fm1,fmid),a2=sim(mid,r,fmid,fm2,fr);
    if(fabs(la-a1-a2)<eps) return a1+a2;
    return calc(l,mid,m1,fl,fmid,fm1,a1)+calc(mid,r,m2,fmid,fr,fm2,a2);
}

int main()
{
    int i,j;
    db l,r;
    cin>>n;
    for(i=1;i<=n;i++)
    {
        scanf("%lf%lf%lf",&tmp[i].x,&tmp[i].y,&tmp[i].r);
    }
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            if(i==j) continue;
            if(len(tmp[i],tmp[j])<tmp[j].r-tmp[i].r) break;
        }
        if(j>n) rd[++nn]=tmp[i];
    }
    n=nn;
    sort(rd+1,rd+n+1);
    for(i=1;i<=n;i++) xl[i]=rd[i].x-rd[i].r,xr[i]=rd[i].x+rd[i].r;
    for(s=1;s<=n;s++)
    {
        l=xl[s],r=xr[s];
        for(t=s+1;t<=n&&r>xl[t];t++) r=max(r,xr[t]);
        t--;
        ans+=calc(l,r,(l+r)/2.0,getf(l),getf(r),getf((l+r)/2.0),1e9);
        s=t;
    }
    printf("%.3f",ans);
}