计算几何收录

题目：
HNOI2016 矿区
SDOI2016 平凡的骰子
JSOI2018 绝地反击
HNOI2004 最佳包裹
SCOI2015 小凸想跑步
APIO2018 选圆圈
lxdAKIMO

HNOI2016 矿区

首先将平面图转成对偶图。将每个点为端点的边按极角排序，然后按逆时针顺序遍历每个区域，算出面积并编号，求出对偶图。
求出对偶图的生成树，一个原图的块在生成树上仍然是联通块，遍历询问的边时判断询问的联通块和该边在树上的关系即可找到对应子树的贡献。
时间复杂度 $O(n{\log }_{2}n+\sum d)$，空间复杂度 $O(n)$

#include<stdio.h>
#include<math.h>
#include<vector>
#include<algorithm>
#include<unordered_map>
using namespace std;
#define R register int
#define L long long
#define I inline
#define N 200001
#define O 400000
#define M 1200000
I L Gcd(L x,L y){
	return y==0?x:Gcd(y,x%y);
}
int PosX[N],PosY[N],f[O],bel[M],rk[M];
bool vis[M];
vector<int>H[O];
L sz2[O],sz1[O],Area[O];
struct Edge{
	int Id,End;
};
I Edge Pair(int y,int d){
	Edge res;
	res.End=y;
	res.Id=d;
	return res;
}
vector<Edge>G[N];
I int GetF(int x){
	if(f[x]==x){
		return x;
	}
	f[x]=GetF(f[x]);
	return f[x];
}
I void DFS(int x,int F){
	sz2[x]=(L)Area[x]*Area[x];
	sz1[x]=Area[x];
	for(int T:H[x]){
		if(T!=F){
			DFS(T,x);
			sz1[x]+=sz1[T];
			sz2[x]+=sz2[T];
		}
	}
}
unordered_map<L,int>Q;
I void Calc(int a,int b,L&res1,L&res2){
	int eid=Q[(L)a*N+b],x,y;
	if(vis[eid]==true){
		x=bel[eid];
		y=bel[eid^1];
		if(sz1[x]>sz1[y]){
			res1-=sz2[y];
			res2-=sz1[y];
		}else{
			res1+=sz2[x];
			res2+=sz1[x];
		}
	}
	
}
int main(){
	int n,m,q,x,y,tot=0,root;
	scanf("%d%d%d",&n,&m,&q);
	for(R i=1;i<=n;i++){
		scanf("%d%d",PosX+i,PosY+i);
	}
	for(R i=0;i!=m;i++){
		scanf("%d%d",&x,&y);
		G[x].push_back(Pair(y,i<<1));
		G[y].push_back(Pair(x,i<<1|1));
		Q[(L)x*N+y]=i<<1;
		Q[(L)y*N+x]=i<<1|1;
	}
	for(R i=1;i<=n;i++){
		vector<pair<double,pair<int,int>>>E;
		for(Edge T:G[i]){
			E.push_back(make_pair(atan2(PosY[T.End]-PosY[i],PosX[T.End]-PosX[i]),make_pair(T.End,T.Id)));
		}
		sort(E.begin(),E.end());
		vector<Edge>().swap(G[i]);
		int r=0;
		for(auto T=E.rbegin();T!=E.rend();T++){
			G[i].push_back(Pair(T->second.first,T->second.second));
			rk[T->second.second]=r;
			r++;
		}
	}
	for(R i=1;i<=n;i++){
		for(Edge T:G[i]){
			if(vis[T.Id]==false){
				L S=0;
				int d1=i,d2=T.End,e=T.Id;
				do{
					vis[e]=true;
					bel[e]=tot;
					S+=PosX[d1]*PosY[d2]-PosY[d1]*PosX[d2];
					d1=d2;
					int r=rk[e^1]+1;
					if(r==G[d2].size()){
						r=0;
					}
					d2=G[d1][r].End;
					e=G[d1][r].Id;
				}while(d1!=i);
				if(S>0){
					Area[tot]=S;
				}else{
					root=tot;
				}
				f[tot]=tot;
				tot++;
			}
		}
	}
	for(R i=0;i!=m;i++){
		x=bel[i<<1];
		y=bel[i<<1|1];
		if(GetF(x)!=GetF(y)){
			f[f[x]]=y;
			H[x].push_back(y);
			H[y].push_back(x);
		}else{
			vis[i<<1]=vis[i<<1|1]=false;
		}
	}
	L ansN=0,ansD;
	Area[root]=1;
	DFS(root,-1);
	for(R i=0;i!=q;i++){
		int d;
		scanf("%d",&d);
		d=(ansN+d)%n+1;
		vector<int>A(d);
		for(R j=0;j!=d;j++){
			scanf("%d",&A[j]);
			A[j]=(ansN+A[j])%n+1;
		}
		ansN=ansD=0;
		Calc(A.back(),A[0],ansN,ansD);
		for(R i=1;i!=d;i++){
			Calc(A[i-1],A[i],ansN,ansD);
		}
		ansD<<=1;
		L g=Gcd(ansN,ansD);
		ansN/=g;
		printf("%lld %lld\n",ansN,ansD/g);
	}
	return 0;
}

平凡的骰子

先要求出骰子的重心。将每个面分成三角形，和原点连成三棱锥，求这些三棱锥各自的重心和体积最后加权平均即可。
统计答案时需要求二面角的大小。二面角用叉乘求出面的法向量求夹角即可。

#include<stdio.h>
#include<math.h>
#include<vector>
using namespace std;
#define R register int
#define D double
#define I inline
#define PIE 3.1415926535897931
struct Vectors{
	D Vx,Vy,Vz;
	I void Read(){
		scanf("%lf%lf%lf",&Vx,&Vy,&Vz);
	}
	I D Size(){
		return sqrt(Vx*Vx+Vy*Vy+Vz*Vz);
	}
	I auto friend operator+(Vectors A,Vectors B){
		A.Vx+=B.Vx;
		A.Vy+=B.Vy;
		A.Vz+=B.Vz;
		return A;
	}
	I auto friend operator-(Vectors A,Vectors B){
		A.Vx-=B.Vx;
		A.Vy-=B.Vy;
		A.Vz-=B.Vz;
		return A;
	}
	I auto friend operator*(Vectors A,D b){
		A.Vx*=b;
		A.Vy*=b;
		A.Vz*=b;
		return A;
	}
	I D friend operator*(Vectors A,Vectors B){
		return A.Vx*B.Vx+A.Vy*B.Vy+A.Vz*B.Vz;
	}
	I auto friend operator^(Vectors A,Vectors B){
		Vectors C;
		C.Vx=A.Vy*B.Vz-A.Vz*B.Vy;
		C.Vy=A.Vz*B.Vx-A.Vx*B.Vz;
		C.Vz=A.Vx*B.Vy-A.Vy*B.Vx;
		return C;
	}
	I D friend operator&(Vectors A,Vectors B){
		return acos(A*B/A.Size()/B.Size());
	}
}p[50];
I auto Pair(D x,D y,D z){
	Vectors res;
	res.Vx=x;
	res.Vy=y;
	res.Vz=z;
	return res;
}
vector<int>A[80];
int main(){
	int n,m,d,x;
	scanf("%d%d",&n,&m);
	for(R i=0;i!=n;i++){
		p[i].Read();
	}
	Vectors G=Pair(0,0,0);
	D TotV=0;
	for(R i=0;i!=m;i++){
		scanf("%d",&d);
		for(R j=0;j!=d;j++){
			scanf("%d",&x);
			A[i].push_back(x-1);
		}
		for(R j=2;j!=d;j++){
			Vectors CurG=p[A[i][0]]+p[A[i][j-1]]+p[A[i][j]];
			D CurV=p[A[i][0]]*(p[A[i][0]]-p[A[i][j]]^p[A[i][j-1]]-p[A[i][0]]);
			G=G+CurG*CurV;
			TotV+=CurV;
		}
	}
	G=G*(.25/TotV);
	for(R i=0;i!=n;i++){
		p[i]=p[i]-G;
	}
	for(R i=0;i!=m;i++){
		D ans=0;
		d=A[i].size();
		Vectors V0=p[A[i][0]];
		for(R j=2;j!=d;j++){
			Vectors V1=p[A[i][j-1]],V2=p[A[i][j]];
			ans+=((V0^V1)&(V2^V1))+((V1^V2)&(V0^V2))+((V2^V0)&(V1^V0))-PIE;
		}
		printf("%.7lf\n",.25*ans/PIE);
	}
	return 0;
}

JSOI2018 绝地反击

根据题意，二分答案。根据答案可以求出每个飞船在圆上可以到达的弧。容易发现，有效的解最多只有 $n$ 个，只要将这 $n$ 个弧的端点逐个尝试即可。确定每个飞船最终的位置之后可以建出二分图判断是否存在完美匹配。由于每个飞船连的点时一个区间，直接断环成链然后用霍尔定理即可。
时间复杂度 $O(n^3\lg R)$，空间复杂度 $O(n)$

#include<stdio.h>
#include<math.h>
#include<vector>
#include<algorithm>
using namespace std;
#define R register int
#define D double
#define I inline
#define EPS 1e-8
#define PIE 3.1415926535897931
I int Min(const int x,const int y){
	return x<y?x:y;
}
int px[200],py[200];
D theL[200],theR[200],w[200];
bool tag[200];
I D Dis(D Ax,D Ay){
	return sqrt(Ax*Ax+Ay*Ay);
}
vector<int>H[400];
I bool CheckRange(D x,D l,D r){
	if(l>r){
		return x>=l||x<=r;
	}
	return l<=x&&x<=r;
}
I bool Solve(int n,D the){
	D del=PIE*2/n;
	int tot=0;
	for(R i=0;i!=n<<1;i++){
		vector<int>().swap(H[i]);
	}
	for(R i=0;i!=n;i++){
		if(tag[i]==true){
			tot++;
		}
		w[i]=the;
		the+=del;
		if(the>=PIE){
			the-=PIE*2;
		}
	}
	sort(w,w+n);
	for(R i=0;i!=n;i++){
		if(tag[i]==false){
			int lf=lower_bound(w,w+n,theL[i])-w,rt;
			if(lf==n){
				lf=0;
			}
			if(CheckRange(w[lf],theL[i],theR[i])==false){
				return false;
			}
			rt=lf;
			int p=rt+1;
			if(p==n){
				p=0;
			}
			while(p!=lf){
				if(CheckRange(w[p],theL[i],theR[i])==false){
					break;
				}
				rt=p;
				p++;
				if(p==n){
					p=0;
				}
			}
			if(lf>rt){
				H[rt+n].push_back(lf);
			}else{
				H[rt].push_back(lf);
				H[rt+n].push_back(lf+n);
			}
		}
	}
	for(R i=0;i<=n;i++){
		int cur=0;
		for(R j=i;j!=i+n;j++){
			for(int T:H[j]){
				if(T>=i){
					cur++;
				}
			}
			if(j==i+n-1){
				cur+=tot;
			}
			if(cur>j-i+1){
				return false;
			}
		}
	}
	return true;
}
I bool Check(int n,int r,D lim){
	vector<D>A;
	for(R i=0;i!=n;i++){
		D l=Dis(px[i],py[i]),cosThe;
		cosThe=(l*l+r*r-lim*lim)/(l*r*2);
		if(cosThe>1){
			return false;
		}
		if(cosThe>-1){
			tag[i]=false;
			D the=atan2(px[i],py[i]),del;
			del=acos(cosThe);
			theL[i]=the-del;
			theR[i]=the+del;
			if(theL[i]<-PIE){
				theL[i]+=PIE*2;
			}
			if(theR[i]>=PIE){
				theR[i]-=PIE*2;
			}
			A.push_back(theL[i]<-PIE?theL[i]+PIE*2:theL[i]);
		}else{
			tag[i]=true;
		}
	}
	if(A.empty()==true){
		A.push_back(0);
	}else{
		sort(A.begin(),A.end());
	}
	D Last=-1e9;
	for(D T:A){
		if(T-Last>EPS){
			if(Solve(n,T)==true){
				return true;
			}
			Last=T;
		}
	}
	return false;
}
int main(){
	int n,r;
	scanf("%d%d",&n,&r);
	for(R i=0;i!=n;i++){
		scanf("%d%d",px+i,py+i);
	}
	D lf=0,rt=300,mid,ans;
	while(lf<rt){
		mid=.5*(lf+rt);
		if(Check(n,r,mid)==true){
			ans=mid;
			rt=mid-EPS;
		}else{
			lf=mid+EPS;
		}
	}
	printf("%.9lf",ans);
	return 0;
}

HNOI2004 最佳包裹

题目要求三维凸包的表面积。先将每个点微扰一下，然后任选四个点出来建成凸包。再逐个加入每个点，加入时如果在某个面之外则删除这个面，最后把删掉的部分补齐即可。存下每个凸包上的面即可实现。
时间复杂度 $O(n^2)$，空间复杂度 $O(n)$

#include<stdio.h>
#include<math.h>
#include<random>
#include<vector>
#include<unordered_set>
#include<time.h>
using namespace std;
#define R register int
#define D double
#define I inline
#define EPS 1e-9
struct Vectors{
	D Vx,Vy,Vz;
	I D Length(){
		return sqrt(Vx*Vx+Vy*Vy+Vz*Vz);
	}
	I auto friend operator-(Vectors A,Vectors B){
		A.Vx-=B.Vx;
		A.Vy-=B.Vy;
		A.Vz-=B.Vz;
		return A;
	}
	I D friend operator*(Vectors A,Vectors B){
		return A.Vx*B.Vx+A.Vy*B.Vy+A.Vz*B.Vz;
	}
	I auto friend operator^(Vectors A,Vectors B){
		Vectors C;
		C.Vx=A.Vy*B.Vz-A.Vz*B.Vy;
		C.Vy=A.Vz*B.Vx-A.Vx*B.Vz;
		C.Vz=A.Vx*B.Vy-A.Vy*B.Vx;
		return C;
	}
}p[100],q[100];
I auto Pair(D x,D y,D z){
	Vectors res;
	res.Vx=x;
	res.Vy=y;
	res.Vz=z;
	return res;
}
struct Faces{
	int VA,VB,VC;
	I void Reverse(){
		int VT=VA;
		VA=VC;
		VC=VT;
	}
	I auto NormalVector(){
		return p[VB]-p[VA]^p[VC]-p[VB];
	}
	I D Area(){
		return NormalVector().Length();
	}
};
I auto MakeFace(int a,int b,int c){
	Faces res;
	res.VA=a;
	res.VB=b;
	res.VC=c;
	return res;
}
I D CalcVolume(Faces F,Vectors V){
	return F.NormalVector()*(V-p[F.VC]);
}
I void Insert(unordered_set<int>&S,int x,int y){
	if(S.count(y<<8|x)!=0){
		S.erase(y<<8|x);
	}else{
		S.insert(x<<8|y);
	}
}
int main(){
	int n;
	scanf("%d",&n);
	mt19937 Rand(time(0));
	for(R i=0;i!=n;i++){
		D x,y,z;
		scanf("%lf%lf%lf",&x,&y,&z);
		q[i]=Pair(x,y,z);
		D a=Rand(),b=Rand();
		p[i]=Pair(cos(a)*cos(b)*EPS+x,sin(a)*cos(b)*EPS+y,sin(b)*EPS+z);
	}
	vector<Faces>F;
	F.push_back(MakeFace(1,2,3));
	F.push_back(MakeFace(0,2,3));
	F.push_back(MakeFace(0,1,3));
	F.push_back(MakeFace(0,1,2));
	for(R i=0;i!=4;i++){
		if(CalcVolume(F[i],p[i])>0){
			F[i].Reverse();
		}
	}
	for(R i=4;i!=n;i++){
		vector<Faces>H;
		unordered_set<int>S;
		for(auto T:F){
			if(CalcVolume(T,p[i])>0){
				Insert(S,T.VA,T.VB);
				Insert(S,T.VB,T.VC);
				Insert(S,T.VC,T.VA);
			}else{
				H.push_back(T);
			}
		}
		F.swap(H);
		for(int T:S){
			int x=T>>8,y=T&255;
			F.push_back(MakeFace(x,y,i));
		}
	}
	for(R i=0;i!=n;i++){
		p[i]=q[i];
	}
	D ans=0;
	for(auto T:F){
		ans+=T.Area();
	}
	printf("%.6lf",.5*ans);
	return 0;
}

SCOI2015 小凸想跑步

根据题意，求出第一条边和其他边的面积平分线，加上第一条边求半平面交面积即可。
时空复杂度 $O(n)$

#include<stdio.h>
#include<math.h>
#include<algorithm>
#include<vector>
using namespace std;
#define R register int
#define L long long
#define D double
#define I inline
#define N 100002
#define EPS 1e-9
struct Point{
	D Vx,Vy;
}p[N];
struct Lines{
	D Vx,Vy,Dx,Dy;
}lim[N],q[N];
I Lines Pair(D x,D y,D a,D b){
	Lines res;
	res.Vx=x;
	res.Vy=y;
	res.Dx=a;
	res.Dy=b;
	return res;
}
I int Check(Lines A,D x,D y){
	x-=A.Vx;
	y-=A.Vy;
	D t=A.Dx*y-x*A.Dy;
	if(t>EPS){
		return 1;
	}
	if(t<-EPS){
		return-1;
	}
	return 0;
}
I auto GetInter(Lines A,Lines B){
	D a=((A.Vy-B.Vy)*B.Dx-(A.Vx-B.Vx)*B.Dy)/(A.Dx*B.Dy-A.Dy*B.Dx);
	Point res;
	res.Vx=A.Vx+a*A.Dx;
	res.Vy=A.Vy+a*A.Dy;
	return res;
}
int px[N],py[N];
int main(){
	int n;
	scanf("%d",&n);
	for(R i=0;i!=n;i++){
		scanf("%d%d",px+i,py+i);
	}
	lim[0]=Pair(px[0],py[0],px[1]-px[0],py[1]-py[0]);
	px[n]=px[0];
	py[n]=py[0];
	L totS=0;
	for(R i=0;i!=n;i++){
		totS+=(L)px[i]*py[i+1]-(L)py[i]*px[i+1];
	}
	for(R i=1;i!=n;i++){
		int ax=px[1]-px[0],ay=py[1]-py[0],bx=px[i+1]-px[i],by=py[i+1]-py[i];
		if((L)ax*by==(L)ay*bx){
			lim[i]=Pair(.5*(px[0]+px[i]),.5*(py[0]+py[i]),bx,by);
		}else{
			Point A=GetInter(lim[0],Pair(px[i],py[i],px[i+1]-px[i],py[i+1]-py[i]));
			lim[i]=Pair(A.Vx,A.Vy,px[0]-px[1]-px[i]+px[i+1],py[0]-py[1]-py[i]+py[i+1]);
			if(Check(lim[i],px[0],py[0])<0||Check(lim[i],px[1],py[1])<0){
				lim[i].Dx=-lim[i].Dx;
				lim[i].Dy=-lim[i].Dy;
			}
		}
	}
	vector<pair<D,int>>C(n);
	for(int i=0;i!=n;i++){
		C[i]=make_pair(atan2(lim[i].Dy,lim[i].Dx),i);
	}
	sort(C.begin(),C.end());
	int head=0,tail=0;
	q[0]=lim[C[0].second];
	for(R i=1;i!=n;i++){
		Lines CurL=lim[C[i].second];
		while(head<tail&&Check(CurL,p[tail].Vx,p[tail].Vy)<1){
			tail--;
		}
		while(head<tail&&Check(CurL,p[head+1].Vx,p[head+1].Vy)<1){
			head++;
		}
		tail++;
		q[tail]=CurL;
		p[tail]=GetInter(q[tail-1],q[tail]);
	}
	while(head<tail&&Check(q[head],p[tail].Vx,p[tail].Vy)<1){
		tail--;
	}
	p[head]=GetInter(q[head],q[tail]);
	D Area=p[tail].Vx*p[head].Vy-p[tail].Vy*p[head].Vx;
	for(R i=head;i!=tail;i++){
		Area+=p[i].Vx*p[i+1].Vy-p[i].Vy*p[i+1].Vx;
	}
	printf("%.4lf",Area/totS);
	return 0;
}

APIO2018 选圆圈

先旋转然后直接建K-D树即可实现删除和查询。
时间复杂度 $O(n\sqrt{n})$，空间复杂度 $O(n)$

#include<stdio.h>
#include<random>
#include<algorithm>
#include<time.h>
using namespace std;
#define R register int
#define L long long
#define D double
#define I inline
#define N 300001
#define INF 2000000001
I void Min(int&x,const int y){
	if(x>y){
		x=y;
	}
}
I void Max(int&x,const int y){
	if(x<y){
		x=y;
	}
}
I bool Check(int lx,int rx,int ly,int ry,int x,int y,int r){
	int dx,dy;
	if(x<lx){
		dx=lx-x;
	}else if(x>rx){
		dx=rx-x;
	}else{
		dx=0;
	}
	if(y<ly){
		dy=ly-y;
	}else if(y>ry){
		dy=y-ry;
	}else{
		dy=0;
	}
	return(L)dx*dx+(L)dy*dy<=(L)r*r;
}
int id[N],ans[N],r[N],pos[N][2];
double p[N][2];
namespace CompareHelper{
	int CurD;
	I bool ComparePos(int x,int y){
		return p[x][CurD^1]<p[y][CurD^1];
	}
	I void Flip(){
		CurD^=1;
	}
	I bool CompareRi(int x,int y){
		if(r[x]==r[y]){
			return x<y;
		}
		return r[x]>r[y];
	}
}
namespace NodeCounter{
	int TotNode;
	I void InitCounter(){
		TotNode=2;
	}
	I void GetNode(int&x){
		x=TotNode;
		TotNode++;
	}
}
struct KDNode{
	int Ls,Rs,Id,Lx,Rx,Ly,Ry;
	I void operator+=(KDNode&A){
		Min(Lx,A.Lx);
		Max(Rx,A.Rx);
		Min(Ly,A.Ly);
		Max(Ry,A.Ry);
	}
}T[N];
I void Update(int p){
	int x=T[p].Id;
	if(x!=0){
		T[p].Lx=pos[x][0]-r[x];
		T[p].Rx=pos[x][0]+r[x];
		T[p].Ly=pos[x][1]-r[x];
		T[p].Ry=pos[x][1]+r[x];
	}else{
		T[p].Lx=T[p].Ly=INF;
		T[p].Rx=T[p].Ry=-INF;
	}
	if(T[p].Ls!=0){
		T[p]+=T[T[p].Ls];
	}
	if(T[p].Rs!=0){
		T[p]+=T[T[p].Rs];
	}
}
I void BuildTree(int p,int lf,int rt){
	int mid=lf+rt>>1;
	nth_element(id+lf,id+mid,id+rt+1,CompareHelper::ComparePos);
	CompareHelper::Flip();
	T[p].Id=id[mid];
	if(mid!=lf){
		NodeCounter::GetNode(T[p].Ls);
		BuildTree(T[p].Ls,lf,mid-1);
	}
	if(mid!=rt){
		NodeCounter::GetNode(T[p].Rs);
		BuildTree(T[p].Rs,mid+1,rt);
	}
	Update(p);
	CompareHelper::Flip();
}
I void Cover(int p,const int Cx,const int Cy,const int Cr,const int Cid){
	if(Check(T[p].Lx,T[p].Rx,T[p].Ly,T[p].Ry,Cx,Cy,Cr)==true){
		if(T[p].Id!=0){
			int dx=Cx-pos[T[p].Id][0],dy=Cy-pos[T[p].Id][1],dr=Cr+r[T[p].Id];
			if((L)dx*dx+(L)dy*dy<=(L)dr*dr){
				ans[T[p].Id]=Cid;
				T[p].Id=0;
			}
		}
		if(T[p].Ls!=0){
			Cover(T[p].Ls,Cx,Cy,Cr,Cid);
			if(T[T[p].Ls].Lx==INF){
				T[p].Ls=0;
			}
		}
		if(T[p].Rs!=0){
			Cover(T[p].Rs,Cx,Cy,Cr,Cid);
			if(T[T[p].Rs].Lx==INF){
				T[p].Rs=0;
			}
		}
	}
}
int main(){
	mt19937 Rand(time(0));
	int n,CX=Rand(),CY=Rand(),x,y;
	scanf("%d",&n);
	double a=Rand(),b;
	b=sin(a);
	a=cos(a);
	for(R i=1;i<=n;i++){
		scanf("%d%d%d",pos[i],pos[i]+1,r+i);
		x=pos[i][0]-CX;
		y=pos[i][1]-CY;
		p[i][0]=a*x-b*y;
		p[i][1]=a*y+b*x;
		id[i]=i;
	}
	NodeCounter::InitCounter();
	BuildTree(1,1,n);
	for(R i=1;i<=n;i++){
		id[i]=i;
	}
	sort(id+1,id+n+1,CompareHelper::CompareRi);
	for(R i=1;i<=n;i++){
		if(ans[id[i]]==0){
			Cover(1,pos[id[i]][0],pos[id[i]][1],r[id[i]],id[i]);
		}
	}
	for(R i=1;i<=n;i++){
		printf("%d ",ans[i]);
	}
	return 0;
}

lxdAKIMO

题目要求椭圆面积异或并，辛普森积分模板。

#include<stdio.h>
#include<vector>
#include<math.h>
#include<algorithm>
using namespace std;
#define R register int
#define D double
#define I inline
I D Square(D x){
	return x*x;
}
I bool CheckEqual(D x,D y){
	x-=y;
	return x>-1e-6&&x<1e-6;
}
struct Oval{
	D VA,VB,VC,VD,VE,VF;
}C[100];
I D Calc(D x,int n){
	vector<D>A;
	for(R i=0;i!=n;i++){
		D a,b,c,del;
		a=C[i].VB;
		b=C[i].VC*x+C[i].VE;
		c=C[i].VA*x*x+C[i].VD*x+C[i].VF;
		del=b*b-a*c*4;
		if(del>0){
			del=sqrt(del);
			A.push_back((del-b)/(a*2));
			A.push_back((-del-b)/(a*2));
		}
	}
	int s=A.size();
	sort(A.begin(),A.end());
	D res=0;
	for(R i=0;i!=s;i+=2){
		res+=A[i+1]-A[i];
	}
	return res;
}
I D Intergral(D lf,D ly,D my,D rt,D ry,int n){
	if(rt-lf<1e-4){
		return(ly+my*4+ry)/6*(rt-lf);
	}
	D vlm=Calc(.25*(lf*3+rt),n),vrm=Calc(.25*(lf+rt*3),n);
	if(CheckEqual(ly+my*4+ry,.5*(ly+ry)+my+(vlm+vrm)*2)==true){
		return(ly+my*4+ry)/6*(rt-lf);
	}
	return Intergral(lf,ly,vlm,.5*(lf+rt),my,n)+Intergral(.5*(lf+rt),my,vrm,rt,ry,n);
}
I D Solve(D l,D r,int n){
	return Intergral(l,Calc(l,n),Calc(.5*(l+r),n),r,Calc(r,n),n);
}
int main(){
	int n;
	scanf("%d",&n);
	vector<pair<D,D>>S;
	for(R i=0;i!=n;i++){
		D a,b,the,x,y,COS,SIN;
		scanf("%lf%lf%lf%lf%lf",&a,&b,&the,&x,&y);
		a=1/(a*a);
		b=1/(b*b);
		COS=cos(the);
		SIN=sin(the);
		C[i].VA=COS*COS*a+SIN*SIN*b;
		C[i].VB=SIN*SIN*a+COS*COS*b;
		C[i].VC=COS*SIN*2*(a-b);
		C[i].VD=((COS*SIN*y-SIN*SIN*x)*b-(COS*SIN*y+COS*COS*x)*a)*2;
		C[i].VE=((COS*SIN*x-COS*COS*y)*b-(COS*SIN*x+SIN*SIN*y)*a)*2;
		C[i].VF=Square(COS*x+SIN*y)*a+Square(SIN*x-COS*y)*b-1;
		a=C[i].VC*C[i].VC-C[i].VA*C[i].VB*4;
		b=C[i].VC*C[i].VE*2-C[i].VB*C[i].VD*4;
		D c=C[i].VE*C[i].VE-C[i].VB*C[i].VF*4,d;
		d=sqrt(b*b-a*c*4);
		D lf=(-d-b)/(a*2),rt=(d-b)/(a*2);
		S.push_back(make_pair(lf<rt?lf:rt,lf>rt?lf:rt));
	}
	sort(S.begin(),S.end());
	D ans=0,lf=S[0].first,rt=S[0].second;
	for(auto T=S.begin()+1;T!=S.end();T++){
		if(T->first>rt){
			ans+=Solve(lf,rt,n);
			lf=T->first;
			rt=T->second;
		}else if(rt<T->second){
			rt=T->second;
		}
	}
	printf("%.2lf",ans+Solve(lf,rt,n));
	return 0;
}