Codechef December Challenge 2018 Before Having Donuts

Before Having Donuts

题目链接

Solution

考虑用一个垂直于$z$轴的平面与截这些甜甜圈，截得的图形是一些圆环的并，设$f(x)$表示过点(0,0,x)垂直于$z$轴的平面截得圆环的面积并，那么答案为$\int f(x)dx$。
我们可以使用自适应辛普森积分去逼近答案。

那现在的问题就是知道$x$，求$f(x)$，也就是求$n$个圆环的面积并。
如果只是圆的面积并那还很好求，直接套个格林公式就可以了。
但现在是圆环，点算？
我们可以用总的面积并减去中间的那些挖空的面积，我们定义顺时针为正方向，逆时针为负方向，我们把每个圆环拆成一个正方向的外圆和一个逆方向的内圆，仔细观察便可以发现，应算正贡献的那块面积外围都是正方向，中间那些应挖空的面积外围都是逆方向，所以只用带上符号，就可以和求圆的面积并一样来做，复杂度$O(n^2)$。

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(in ti=j;i>=l;--i)

using namespace std;
typedef long long ll;
typedef double db;

const ll N=160;
const db pi=acos(-1),LE=-1000,RI=800;
const db eps=1e-9;

db x[N],y[N],z[N],R[N],r[N];
db _R[N];

db ang[N];
int xs[N];
int re[N];

int n,m,t,s,k;

inline void join(db R,db x,db y,db r,int val)
{
	db dis=sqrt(x*x+y*y);
	if(dis>=R+r||dis+r<=R)return;
	if(R+dis<=r){
		ang[++m]=0; xs[m]=val;
		ang[++m]=2*pi; xs[m]=-val;
		return;
	}
	db op=acos((R*R+dis*dis-r*r)/(2*R*dis));
	db an=atan2(y,x);
	db lef=an-op,rig=an+op;
	if(rig<0){
		rig=rig+2*pi;
		if(lef<0&&lef+2*pi<=rig)lef+=2*pi;
	}
	if(lef<0){
		ang[++m]=0; xs[m]=val;
		ang[++m]=rig; xs[m]=-val;
		ang[++m]=lef+2*pi; xs[m]=val;
		ang[++m]=2*pi; xs[m]=-val;
	} else {
		ang[++m]=lef; xs[m]=val;
		ang[++m]=rig; xs[m]=-val;
	}
}

inline bool kmp(int a,int b)
{return ang[a]<ang[b];}

inline db board(int bh,db _x,db _y,db p)
{
	m=0;
	fo(i,1,n)if(_R[i]!=-1&&i!=bh)
	join(p,x[i]-_x,y[i]-_y,R[i]+_R[i],1),join(p,x[i]-_x,y[i]-_y,R[i]-_R[i],-1);
	ang[++m]=2*pi; xs[m]=0;
	fo(i,1,m)re[i]=i;
	db las=0;
	int qz=0;
	sort(re+1,re+m+1,kmp);
	db a1=0,a2=0;
	fo(i,1,m){
		int k=re[i];
		if(!qz){
			a1=a1+ang[k]-las;
			a2=a2+(sin(ang[k])-sin(las))*_x+(cos(las)-cos(ang[k]))*_y;
		}
		las=ang[k],qz+=xs[k];
	}
	db ans=a1*p*p+a2*p;
	return ans;
}

inline db cal(db h)
{
	fo(i,1,n)if(abs(h-z[i])<=r[i])_R[i]=sqrt(r[i]*r[i]-(h-z[i])*(h-z[i]));
	else _R[i]=-1;
	db ans=0;
	fo(i,1,n)if(_R[i]!=-1)
	ans=ans+board(i,x[i],y[i],R[i]+_R[i])-board(i,x[i],y[i],R[i]-_R[i]);
	return ans/2;
}

inline db get(db le,db ri)
{return (ri-le)*(cal(le)+cal(ri)+cal((le+ri)*0.5)*4.000)/6.00;}

inline db ccc(db le,db ri,db gg)
{
	db mid=(le+ri)*0.5;
	db fr=get(le,mid),ba=get(mid,ri);
	if(abs(fr+ba-gg)<=eps&&ri-le<=0.3)return gg;
	return ccc(le,mid,fr)+ccc(mid,ri,ba);
}

int main()
{
	scanf("%d",&n);
	fo(i,1,n)scanf("%lf%lf%lf%lf%lf",&x[i],&y[i],&z[i],&R[i],&r[i]);
	db gg=get(LE,RI);
	db ans=ccc(LE,RI,gg);
	ans=ans;
	printf("%.5lf",ans);
}