圆的面积并

Posted 2020-12-17 bestlxm

以前以为圆的面积并只能辛普森积分。今天才知道原来可以(O(n^2log n))做。
大致思路是对每个圆求出不与其他圆重合的角度区域，加上这个角度对应的弓形面积，并且加上这段区域边界两个点的叉积的(frac{1}{2})。
仔细画一画图发现这恰好是对的。非常神奇qwq
代码如下

double sqr(double x){return x*x;}
struct vec{
	double x,y;
	vec operator+ (const vec b) const{return (vec){x+b.x,y+b.y};}
	vec operator- (const vec b) const{return (vec){x-b.x,y-b.y};}
	vec operator* (const double b) const{return (vec){x*b,y*b};}
	double operator* (const vec b) const{return x*b.x+y*b.y;}
	double operator^ (const vec b) const{return x*b.y-y*b.x;}
	double norm(){return sqrt(sqr(x)+sqr(y));}
	double norm2(){return sqr(x)+sqr(y);}
	double ang(){return atan2(y,x);}
}P[M];
vec go(vec x,double k,double r){return (vec){x.x+r*cos(k),x.y+r*sin(k)};}
double calc(int key){
	if (R[key]<1e-9) return 0;double ans=0;int sum=0;
	static std::pair <double,int> Tag[M];int cnt=0;
	for (register int i=1;i<=n;i++){
		if (i==key||R[i]<1e-9) continue;
		double len=(P[i]-P[key]).norm();
		if (R[i]<R[key]&&len<=R[key]-R[i]) continue;
		if (R[i]>R[key]&&len<=R[i]-R[key]) return 0;
		if (R[i]+R[key]<=len) continue;
		double mid=(P[i]-P[key]).ang();
		double ang=acos((sqr(R[key])+sqr(len)-sqr(R[i]))/(2*len*R[key]));
		double ml=mid-ang,mr=mid+ang;if (ml<-pi) ml+=2*pi;if (mr>pi) mr-=2*pi;
		Tag[++cnt]=std::make_pair(ml,1),Tag[++cnt]=std::make_pair(mr,-1);
		if (ml>mr) ++sum;
	}
	Tag[++cnt]=std::make_pair(-pi,0),Tag[++cnt]=std::make_pair(pi,1);
	std::sort(Tag+1,Tag+1+cnt);
	for (register int i=1,j=1;i<=cnt;i=j){
		while (j<=cnt&&Tag[j].first-Tag[i].first<1e-8) sum+=Tag[j++].second;
		if (!sum){
			double ang=Tag[j].first-Tag[i].first;
			ans+=sqr(R[key])*(ang-sin(ang))+(go(P[key],Tag[i].first,R[key])^go(P[key],Tag[j].first,R[key]));
		}
	}
	return ans;
}
double Sum(){double ans=0;for (register int i=1;i<=n;i++) ans+=calc(i);return ans/2;}

这个东西再套个辛普森积分就可以求球的面积并了。