SPOJ 8073 CIRU 自适应Simpson积分求圆的面积并 解题报告

CIRU - The area of the union of circles

You are given N circles and expected to calculate the area of the union of the circles !

Input

The first line is one integer n indicates the number of the circles. (1 <= n <= 1000)
Then follows n lines every line has three integers
Xi Yi Ri
indicates the coordinate of the center of the circle, and the radius. (|Xi|. |Yi| <= 1000, Ri <= 1000)
Note that in this problem Ri may be 0 and it just means one point !

Output

The total area that these N circles with 3 digits after decimal point

Example

Input:
3
0 0 1
0 0 1
100 100 1
Output:
6.283
Submit solution!

【解题报告】
原来以为这个东西很高深，最后发现也不过如此。。（虽然还是记的公式，不够要记得就真的只有一句）
这是一道求圆面积并的题，我们采用自适应Simpson积分来求它的面积（也可以用纯几何的方法，不过个人感觉那种方法不灵活,推起来复杂）。
首先，这些圆应该是长成一坨一坨的（有一些圆互相部分重叠变成了一坨），那么我们的自适应积分就是用来处理这连续的一坨的，所以先去掉完全包含的情况后（即prework），我们用f(x)来表示一根穿过了横坐标为x的扫描线，被圆形覆盖的长度，那么f(x)的图像是在一段一段的区间中是连续的，并且函数的曲线肯定是在x轴的上方。
而且显然最终的圆并面积就是f(x)的图像与x轴圈起来的面积，我们可以用Simpson积分来求这个图像的面积，那么到底什么是Simpson积分呢？它的大概原理就是我们取这条f(x)上的三个点，用这过三个点的二次函数来代替原有曲线，便可以方便快速地求出它的面积，精度最高。
那么什么是自适应呢？简单地说就是在当前区间精度不优的话就二分递归到左右两个区间继续求，如果直接拟合的精度足够的情况下，就不会继续递归到左右两半区间再做Simpson积分，提高了算法的效率。

代码如下：

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define N 10010
#define eps 1e-6
#define inf 0x3f3f3f3f

struct Point
{  
    double x,y;         
    Point(double _=.0,double __=.0):x(_),y(__){}  
    friend bool operator <(const Point &a,const Point &b)
    {return a.x<b.x||(a.x==b.y&&a.y<b.y);}
    Point operator +(const Point &a)const 
    {return Point(x+a.x,y+a.y);}  
    Point operator -(const Point &a)const 
    {return Point(x-a.x,y-a.y);}  
    Point operator *(double a)const 
    {return Point(x*a,y*a);}
};
struct Circle
{
    Point o;
    double r;
    Circle(Point _o=.0,double _r=.0):o(_o),r(_r){}
}circles[N],tmp[N];

bool cmp2(Circle a,Circle b){return a.r>b.r;}
bool cmp(Circle a,Circle b)
{return(fabs(a.o.x-a.r-b.o.x+b.r)<=eps)?a.o.x+a.r<b.o.x+b.r:a.o.x-a.r<b.o.x-b.r;}
bool check(Circle a,Circle b) 
{return (a.o.x-b.o.x)*(a.o.x-b.o.x)+(a.o.y-b.o.y)*(a.o.y-b.o.y)<=(a.r-b.r)*(a.r-b.r);}
void prework(int n)
{
    sort(circles+1,circles+n+1,cmp2);
    int cnt=0;
    for(int i=1;i<=n;++i)
    {
        bool flag=true;
        for(int j=1;j<=cnt;++j)
        if(check(tmp[j],circles[i]))
        {
            flag=false;
            break;
        }
        if(flag) tmp[++cnt]=circles[i];
    }
    n=cnt;
    memcpy(circles,tmp,sizeof(tmp));
}

pair<double,double>intervals[N];
int tot=0,st,ed;

void getCircleIntersec(Circle a,double x)
{
    double len=sqrt(a.r*a.r-(x-a.o.x)*(x-a.o.x));
    intervals[++tot]=make_pair(a.o.y-len,a.o.y+len);
}
double f(double x) 
{
    tot=0;
    for(int i=st;i<=ed;i++)
    if(x<circles[i].o.x+circles[i].r&&x>circles[i].o.x-circles[i].r)
        getCircleIntersec(circles[i],x);
    sort(intervals+1,intervals+tot+1);
    double ans=0,start=-inf,end=-inf; 
    for(int i=1;i<=tot;++i)
    {
        if(intervals[i].first>=end)
        {
            ans+=end-start;
            start=intervals[i].first;
            end=intervals[i].second;
        }
        else end=max(end,intervals[i].second);
    }
    ans+=end-start;
    return ans;
}
double calc(double len,double fL,double fM,double fR) 
{return (fL+4*fM+fR)*len/6;}
double Simpson(double L,double M,double R,double fL,double fM,double fR,double sqr) 
{
    double M1=(L+M)/2,M2=(M+R)/2;
    double fM1=f(M1),fM2=f(M2);
    double g1=calc(M-L,fL,fM1,fM),g2=calc(R-M,fM,fM2,fR);
    if(fabs(sqr-g1-g2)<=eps) return g1+g2;
    return Simpson(L,M1,M,fL,fM1,fM,g1)+Simpson(M,M2,R,fM,fM2,fR,g2);
}

int main()
{
    double ans=0;int n;
    scanf("%d",&n);
    for(int i=1;i<=n;++i) scanf("%lf%lf%lf",&circles[i].o.x,&circles[i].o.y,&circles[i].r);
    prework(n);
    sort(circles+1,circles+n+1,cmp);
    for(int i=1,j;i<=n;++i)
    {
        double L=circles[i].o.x-circles[i].r,R=circles[i].o.x+circles[i].r;
        for(j=i+1;j<=n;j++)
        {
            if(circles[j].o.x-circles[j].r>R) break;
            else R=max(R,circles[j].o.x+circles[j].r); 
        }
        double M=(L+R)/2;
        st=i,ed=j-1;
        i=j-1;
        double fL=f(L),fM=f(M),fR=f(R);
        ans+=Simpson(L,M,R,fL,fM,fR,calc(R-L,fL,fM,fR));
    }
    printf("%.3lf\n",ans);
    return 0;
}