[SPOJ CIRU]The area of the union of circles(自适应Simpson积分求圆并面积)

题目链接

http://acm.hust.edu.cn/vjudge/problem/viewProblem.action?id=16754

题目大意

求圆的面积并

思路

首先膜拜一发福大核武的神做法
http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6
这种做法在applepi大爷的计算几何讲义中也有提到过，非常神奇优美，但是很不好写，毕竟不是所有人都像核武大爷那样神犇，没有很强的可推广性。

再来膜拜一发七爷的自适应Simpson积分题解
http://hi.baidu.com/sevenkplus/item/8cf47779e1e83b376cc37caa

自适应Simpson积分来求圆并的做法，就是用$f(x)$来表示一根穿过了横坐标为$x$的扫描线，被圆形覆盖的长度，那么$f(x)$的图像是在一段一段的区间中是连续的，并且函数的曲线肯定是在$x$轴的上方。而且显然最终的圆并面积就是$f(x)$的图像与$x$轴圈起来的面积，我们可以用数值积分来求这个图像的面积，数值积分分三种：矩形、梯形、Simpson，其中Simpson积分的做法就是把函数图像用一根二次函数曲线来拟合，精度最高。
而自适应Simpson积分，简单地说就是在当前区间直接拟合的精度足够的情况下，就不会继续递归到左右两半区间再做Simpson积分，提高了算法的效率。

看起来这个玩意就是乱搞嘛是不是= =，似乎看上去很不靠谱，但是实际上非常好用，除非一些特殊的反例(http://wapapp.baidu.com/shuxk/item/bd1210284b5d2e869c63d1b7)，不光是圆的面积并，只要你会求出某种图形盖在扫描线上的区间，几乎任何玩意的面积并都是可以求的，实为杀人越货之利器23333。

那么我们就是要把所有圆分成一坨一坨圆，每一坨圆是连在一块的，对于每一坨圆，它们的$f(x)$图像肯定是连续的，就能用自适应Simpson积分把这一坨圆的面积求出来，累加答案即可。
当然还要注意怎么求一根扫描线被圆覆盖的长度，这和之前我做的下落的圆盘那题差不多，也是个比较简单的区间贪心问题，而找一坨一坨圆也是一样的贪心，只不过是在$x$轴上做贪心而已，怎么做相信大家都会，渣渣我就不说了。

嘴巴上讲讲做法当然很简单，但是我还是花了一个白天的时间才A掉这题，因此还要注意一些细节：
1、EPS开1e-6足矣，1e-7也可以，不过速度慢了4倍
2、此题很卡精度和细节，判$<$EPS神马的必须取等号，比如说写成 < <script id="MathJax-Element-1121" type="math/tex"><</script>EPS就会WA，写成 <= <script id="MathJax-Element-1122" type="math/tex"><=</script>EPS就没问题，不过在1e-7精度下没问题，坑爹啊
3、找区间并的长度的贪心很简单，但是还是要注意，不能掉以轻心（看我注释中的感叹号就是我WA的地方）
4、要做一个预处理，删掉所有被其他大圆完全覆盖的小圆，否则会WA，不要问我为什么，我也不知道，亲测WA，扒了别人代码才发现问题的。

代码

//自适应Simpson积分
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <cmath>

#define MAXN 10000
#define EPS 1e-6

using namespace std;

int n;

struct Point //点
{
    double x,y;
    Point(){}
    Point(double _x,double _y):x(_x),y(_y){}
    bool operator<(const Point &b)const { return x==b.x?y<b.y:x<b.x; }
};

struct Circle //圆
{
    Point o;
    double r;
    Circle(){}
    Circle(Point _o,double _r):o(_o),r(_r){}
}circles[MAXN],tmp[MAXN];

bool cmp2(Circle a,Circle b)
{
    return a.r>b.r;
}

bool cmp(Circle a,Circle b)
{
    if(fabs(a.o.x-a.r-b.o.x+b.r)<=EPS) return a.o.x+a.r<b.o.x+b.r;
    return a.o.x-a.r<b.o.x-b.r;
}

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

bool check(Circle a,Circle b) //检查b是否被套在a里面
{
    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()
{
    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;
    for(int i=1;i<=n;i++)
        circles[i]=tmp[i];
}

void getCircleIntersec(Circle a,double x) //求x=x这条平行于y轴的直线穿过圆a的长度
{
    intervals[++tot]=make_pair(a.o.y-sqrt(a.r*a.r-(x-a.o.x)*(x-a.o.x)),a.o.y+sqrt(a.r*a.r-(x-a.o.x)*(x-a.o.x)));
}

double f(double x) //求扫描线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=-1e12,end=-1e12; //!!!!!!!!!!!
    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) //求长度为len的[L,R]区间，中点为M的Simpson近似面积
{
    return (fL+4*fM+fR)*len/6;
}

double Simpson(double L,double M,double R,double fL,double fM,double fR,double sqr) //Simpson积分求区间[L,R]的面积并,f(L)=L,f(R)=R,f(M)=M,把[L,R]当成整体来拟合得到的面积是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) //把当前区间分成2半再拟合得到的答案差别很小，就不再递归下去了
        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;
    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();
    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;
}