[SPOJ-CIRU]The area of the union of circles/[BZOJ2178]圓的面積並

題目大意：

求\(n(n\le1000)\)個圓的面積並。

思路：

對於一個\(x\)，我們可以用線段覆蓋的方法求出被圓覆蓋的長度。用\(f(x)\)表示橫座標為\(x\)時覆蓋的長度，則我們可以對\(f(x)\)積分來得到答案。注意面積不連續的部分要分開求。

原始碼：

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<algorithm>
inline int getint() {
    register char ch;
    register bool neg=false;
    while(!isdigit(ch=getchar())) neg|=ch=='-';
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return neg?-x:x;
}
const int N=1e3+1;
const double eps=1e-4;
inline double sqr(const double &x) {
    return x*x;
}
int n;
struct Circle {
    double x,y,r;
};
Circle c[N];
struct Segment {
    double l,r;
    bool operator < (const Segment &rhs) const {
        return r<rhs.r;
    }
};
Segment t[N];
bool mark[N];
std::vector<std::pair<double,int> > s;
inline double F(const double &x) {
    for(register int i=1;i<=n;i++) {
        if(fabs(c[i].x-x)>c[i].r) continue;
        const double d=sqrt(sqr(c[i].r)-sqr(c[i].x-x));
        s.push_back(std::make_pair(c[i].y-d,1));
        s.push_back(std::make_pair(c[i].y+d,-1));
    }
    std::sort(s.begin(),s.end());
    int tmp=0;
    double st,ans=0;
    for(register unsigned i=0;i<s.size();i++) {
        if(tmp==0) st=s[i].first;
        tmp+=s[i].second;
        if(tmp==0) ans+=s[i].first-st;
    }
    s.clear();
    return ans;
}
inline double simpson(const double &a,const double &b) {
    const double c=(a+b)/2;
    return (F(a)+4*F(c)+F(b))*(b-a)/6;
}
inline double asr(const double &a,const double &b,const double &eps,const double &s) {
    const double &c=(a+b)/2,ls=simpson(a,c),rs=simpson(c,b);
    if(fabs(ls+rs-s)<=15*eps) return ls+rs+(ls+rs-s)/15;
    return asr(a,c,eps/2,ls)+asr(c,b,eps/2,rs);
}
int main() {
    n=getint();
    for(register int i=1;i<=n;i++) {
        c[i].x=getint();
        c[i].y=getint();
        c[i].r=getint();
    }
    for(register int i=1;i<=n;i++) {
        for(register int j=1;j<=n;j++) {
            if(c[i].r>c[j].r||i==j) continue;
            if(sqr(c[i].x-c[j].x)+sqr(c[i].y-c[j].y)<sqr(c[i].r-c[j].r)) {
                mark[i]=true;
            }
        }
    }
    int tot=0;
    for(register int i=1;i<=n;i++) {
        if(!mark[i]) c[++tot]=c[i];
    }
    n=tot;
    for(register int i=1;i<=n;i++) {
        t[i]=(Segment){c[i].x-c[i].r,c[i].x+c[i].r};
    }
    std::fill(&mark[1],&mark[n]+1,0);
    for(register int i=1;i<=n;i++) {
        for(register int j=1;j<=n;j++) {
            if(i==j) continue;
            if(t[j].l<t[i].l&&t[i].r<t[j].r) mark[i]=true;
        }
    }
    tot=0;
    for(register int i=1;i<=n;i++) {
        if(!mark[i]) t[++tot]=t[i];
    }
    std::sort(&t[1],&t[tot]+1);
    double st=t[1].l,en=t[1].r,ans=0;
    for(register int i=2;i<=tot;i++) {
        if(t[i].l>en) {
            ans+=asr(st,en,eps,simpson(st,en));
            st=t[i].l;
        }
        en=t[i].r;
    }
    ans+=asr(st,en,eps,simpson(st,en));
    printf("%.3f\n",ans);
    return 0;
}
 

但是BZOJ上比較卡精度，而eps開太小又會T，解決方法是將Lyness優化去掉，將沒有用的圓去掉。

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<algorithm>
inline int getint() {
    register char ch;
    register bool neg=false;
    while(!isdigit(ch=getchar())) neg|=ch=='-';
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return neg?-x:x;
}
const int N=1e3+1;
const double eps=1e-13;
inline double sqr(const double &x) {
    return x*x;
}
int n;
struct Circle {
    double x,y,r;
};
Circle c[N];
bool mark[N];
std::vector<int> cir;
std::vector<std::pair<double,int> > s;
inline double F(const double &x) {
    for(register unsigned j=0;j<cir.size();j++) {
        const int &i=cir[j];
        if(std::abs(c[i].x-x)>c[i].r) continue;
        const double d=sqrt(sqr(c[i].r)-sqr(c[i].x-x));
        s.push_back(std::make_pair(c[i].y-d,1));
        s.push_back(std::make_pair(c[i].y+d,-1));
    }
    std::sort(s.begin(),s.end());
    double st,ans=0;
    for(register unsigned i=0,tmp=0;i<s.size();i++) {
        if(tmp==0) st=s[i].first;
        tmp+=s[i].second;
        if(tmp==0) ans+=s[i].first-st;
    }
    s.clear();
    return ans;
}
inline double simpson(const double &fl,const double &fm,const double &fr,const double &len) {
    return (fl+4*fm+fr)*len/6;
}
inline double asr(const double &a,const double &b,const double &fl,const double &fm,const double &fr) {
    const double c=(a+b)/2;
    const double flm=F((a+c)/2),frm=F((c+b)/2);
    const double ls=simpson(fl,flm,fm,c-a),rs=simpson(fm,frm,fr,b-c),s=simpson(fl,fm,fr,b-a);
    if(std::abs(ls+rs-s)<eps) return ls+rs;
    return asr(a,c,fl,flm,fm)+asr(c,b,fm,frm,fr);
}
struct Node {
    double p;
    int v,id;
    bool operator < (const Node &rhs) const {
        return p<rhs.p;
    }
};
std::vector<Node> t;
inline double dist(const Circle &a,const Circle &b) {
    return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
int main() {
    n=getint();
    for(register int i=1;i<=n;i++) {
        c[i].x=getint();
        c[i].y=getint();
        c[i].r=getint();
    }
    for(register int i=1;i<=n;i++) {
        for(register int j=1;j<=n;j++) {
            if(i==j) continue;
            if(dist(c[i],c[j])<c[j].r-c[i].r) {
                mark[i]=true;
                break;
            }
        }
    }
    int tot=0;
    for(register int i=1;i<=n;i++) {
        if(!mark[i]) c[++tot]=c[i];
    }
    n=tot;
    for(register int i=1;i<=n;i++) {
        t.push_back((Node){c[i].x-c[i].r,1,i});
        t.push_back((Node){c[i].x+c[i].r,-1,i});
    }
    std::sort(t.begin(),t.end());
    int beg;
    double ans=0;
    for(register unsigned i=0,tmp=0;i<t.size();i++) {
        if(tmp==0) beg=i;
        tmp+=t[i].v;
        if(tmp==0) {
            for(register int j=beg;j<=(int)i;j++) {
                cir.push_back(t[j].id);
            }
            std::sort(cir.begin(),cir.end());
            cir.resize(std::unique(cir.begin(),cir.end())-cir.begin());
            const double &st=t[beg].p,&en=t[i].p,mid=(st+en)/2;
            ans+=asr(st,en,F(st),F(mid),F(en));
            cir.clear();
        }
    }
    printf("%.3f\n",ans);
    return 0;
}