[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神马的必须取等号，比如说写成<EPS就会WA，写成<=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-6using 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;}