[SPOJ-CIRU]The area of the union of circles/[BZOJ2178]圆的面积并

[SPOJ-CIRU]The area of the union of circles/[BZOJ2178]圆的面积并

题目大意：

求(n(nle1000))个圆的面积并。

思路：

对于一个(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
",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
",ans);
	return 0;
}