计算几何 / Geometry（ACM-Template 2.0）

计算几何

struct of 向量

  • rotate()返回逆时针旋转后的点，left()返回朝左的单位向量
  • trans()返回p沿a,b拉伸的结果，arctrans()返回p在坐标系<a,b>中的坐标
  • 常量式写法，不要另加变量，需要加变量就再搞个struct
  • 直线类在半面交里，其中包含线段交点
struct vec{
	lf x,y; vec(){} vec(lf x,lf y):x(x),y(y){}
	vec operator-(const vec &b){return vec(x-b.x,y-b.y);}
	vec operator+(const vec &b){return vec(x+b.x,y+b.y);}
	vec operator*(lf k){return vec(k*x,k*y);}
	lf len(){return hypot(x,y);}
	lf sqr(){return x*x+y*y;}
	vec trunc(lf k=1){return *this*(k/len());}
	vec rotate(double th){lf c=cos(th),s=sin(th); return vec(x*c-y*s,x*s+y*c);}
	vec left(){return vec(-y,x).trunc();}
	lf theta(){return atan2(y,x);}
	friend lf cross(vec a,vec b){return a.x*b.y-a.y*b.x;};
	friend lf cross(vec a,vec b,vec c){return cross(a-c,b-c);}
	friend lf dot(vec a,vec b){return a.x*b.x+a.y*b.y;}
	friend vec trans(vec p,vec a,vec b){
		swap(a.y,b.x);
		return vec(dot(a,p),dot(b,p));
	}
	friend vec arctrans(vec p,vec a,vec b){
		lf t=cross(a,b);
		return vec(-cross(b,p)/t,cross(a,p)/t);
	}
	void output(){printf("%.12f %.12f\n",x,y);}
}a[N];
vec projection(vec v,vec a,vec b){ // v 在 line(a,b) 上的投影
	vec d=b-a;
	return a+d*(dot(v-a,d)/d.sqr());
}
  • 整数向量
struct vec{
	int x,y; vec(){} vec(int x,int y):x(x),y(y){}
	vec operator-(const vec &b){return vec(x-b.x,y-b.y);}
	vec operator+(const vec &b){return vec(x+b.x,y+b.y);}
	void operator+=(const vec &b){x+=b.x,y+=b.y;}
	void operator-=(const vec &b){x-=b.x,y-=b.y;}
	vec operator*(lf k){return vec(k*x,k*y);}
	bool operator==(vec b)const{return x==b.x && y==b.y;}
	int sqr(){return x*x+y*y;}
	void output(){printf("%lld %lld\n",x,y);}
}a[N]; const vec dn[]={{1,0},{0,1},{-1,0},{0,-1},{1,1},{1,-1},{-1,1},{-1,-1}};

struct of 线段

struct line{
	vec p1,p2; lf th;
	line(){}
	line(vec p1,vec p2):p1(p1),p2(p2){
		th=(p2-p1).theta();
	}
	bool contain(vec v){ // contains the point
		return cross(v,p2,p1)<=eps;
	}
	vec PI(line b){ // point of intersection
		lf t1=cross(p1,b.p2,b.p1);
		lf t2=cross(p2,b.p2,b.p1);
		return vec((t1*p2.x-t2*p1.x)/(t1-t2),(t1*p2.y-t2*p1.y)/(t1-t2));
	}
};

struct of 圆

struct cir{
	vec v; lf r;
	void PI(vec a,vec b,vec &A,vec &B){ // PI with line(a,b)
		vec H=projection(v,a,b);
		vec D=(a-b).trunc(sqrt(r*r-(v-H).sqr()));
		A=H+D; B=H-D;
	}
	void PI(cir b,vec &A,vec &B){ // PI with circle
		vec d=b.v-v;
		lf dis=abs(b.r*b.r-r*r-d.sqr())/(2*d.len());
		vec H=v+d.trunc(dis);
		vec D=d.left().trunc(sqrt(r*r-dis*dis));
		A=H+D; B=H-D;
	}
};

平面几何基本操作

判断两条线段是否相交

  • 快速排斥实验：判断线段所在矩形是否相交（用来减小常数，可省略）
  • 跨立实验：任一线段的两端点在另一线段的两侧
bool judge(vec a,vec b,vec c,vec d){ // 线段ab和线段cd
	#define SJ(x) max(a.x,b.x)<min(c.x,d.x)\
	|| max(c.x,d.x)<min(a.x,b.x)
	if(SJ(x) || SJ(y))return 0;
	#define SJ2(a,b,c,d) cross(a-b,a-c)*cross(a-b,a-d)<=0
	return SJ2(a,b,c,d) && SJ2(c,d,a,b);
}

点是否在线段上

bool onseg(vec p,vec a,vec b){
	return (a.x-p.x)*(b.x-p.x)<eps
	&& (a.y-p.y)*(b.y-p.y)<eps
	&& abs(cross(a-b,a-p))<eps;
}

多边形面积

lf area(vec a[],int n){
	lf ans=0;
	repeat(i,0,n)
		ans+=cross(a[i],a[(i+1)%n]);
	return abs(ans/2);
}

多边形的面积质心

vec centre(vec a[],int n){
	lf S=0; vec v=vec();
	repeat(i,0,n){
		vec &v1=a[i],&v2=a[(i+1)%n];
		lf s=cross(v1,v2);
		S+=s; v=v+(v1+v2)*s;
	}
	return v*(1/(3*S));
}

二维凸包

  • 求上凸包，按坐标 (x, y) 字典升序排序，从小到大加入栈，如果出现凹多边形情况则出栈。下凸包反着来
  • \(O(n\log n)\)，排序是瓶颈
vector<vec> st;
void push(vec &v,int b){
	while((int)st.size()>b
	&& cross(*++st.rbegin(),st.back(),v)<=0) // 会得到逆时针的凸包
		st.pop_back();
	st.push_back(v);
}
void convex(vec a[],int n){
	st.clear();
	sort(a,a+n,[](vec a,vec b){
		return make_pair(a.x,a.y)<make_pair(b.x,b.y);
	});
	repeat(i,0,n)push(a[i],1);
	int b=st.size();
	repeat_back(i,1,n-1)push(a[i],b); // repeat_back自动变成上凸包
}

<补充> 动态凸包

  • 支持添加点、询问点是否在凸包内，\(O(\log n)\)
const lf eps=1e-7;
multiset<pair<lf,vec>> st; vec c;
typedef multiset<pair<lf,vec>>::iterator ptr;
void push(vec v){
	st.insert({v.theta(),v});
}
void dec(ptr &p){if(p==st.begin())p=st.end(); --p;}
void inc(ptr &p){++p; if(p==st.end())p=st.begin();}
ptr find(vec v){
	auto p=st.lower_bound({v.theta(),v}); dec(p);
	return p;
}
bool out(vec v){ // whether out of the convex
	v=v-c;
	auto l=find(v),r=l; inc(r);
	return cross(l->se,r->se,v)<-eps;
}
void init(vec v1,vec v2,vec v3){
	st.clear();
	c=(v1+v2+v3)*(1.0/3);
	push(v1-c); push(v2-c); push(v3-c);
}
void add(vec v){ // add a point to convex
	if(!out(v))return;
	v=v-c;
	auto l=find(v),r=l; inc(r);
	auto l2=l; dec(l2);
	while(cross(l->se,l2->se,v)>0){
		dec(l),dec(l2);
	}
	auto r2=r; inc(r2);
	while(cross(r->se,r2->se,v)<0){
		inc(r),inc(r2);
	}
	inc(l);
	if(l<=r)st.erase(l,r);
	else st.erase(l,st.end()),st.erase(st.begin(),r);
	st.insert({v.theta(),v});
}

旋转卡壳

  • 每次找到凸包每条边的最远点，基于二维凸包，\(O(n\log n)\)
lf calipers(vec a[],int n){
	convex(a,n); // 凸包算法
	repeat(i,0,st.size())a[i]=st[i]; n=st.size();
	lf ans=0; int p=1; a[n]=a[0];
	repeat(i,0,n){
		while(cross(a[p],a[i],a[i+1])<cross(a[p+1],a[i],a[i+1])) // 必须逆时针凸包
			p=(p+1)%n;
		ans=max(ans,(a[p]-a[i]).len());
		ans=max(ans,(a[p+1]-a[i]).len()); // 这里求了直径
	}
	return ans;
}

最大空矩形 using 扫描法

  • 在范围 (0, 0) 到 (l, w) 内求面积最大的不覆盖任何点的矩形面积，\(O(n^2)\)，n 是点数
  • 如果是 lf 就把 vec 结构体内部、ans、u和 d 的类型改一下
struct vec{
	int x,y; // 可能是lf
	vec(int x,int y):x(x),y(y){}
};
vector<vec> a; // 存放点
int l,w;
int ans=0;
void work(int i){
	int u=w,d=0;
	repeat(k,i+1,a.size())
	if(a[k].y>d && a[k].y<u){
		ans=max(ans,(a[k].x-a[i].x)*(u-d)); // 更新ans
		if(a[k].y==a[i].y)return; // 可行性剪枝
		(a[k].y>a[i].y?u:d)=a[k].y; // 更新u和d
		if((l-a[i].x)*(u-d)<=ans)return; // 最优性剪枝
	}
	ans=max(ans,(l-a[i].x)*(u-d)); // 撞墙更新ans
}
int query(){
	a.push_back(vec(0,0));
	a.push_back(vec(l,w)); // 加两个点方便处理
	// 小矩形的左边靠着顶点的情况
	sort(a.begin(),a.end(),[](vec a,vec b){return a.x<b.x;});
	repeat(i,0,a.size())
		work(i);
	// 小矩形的右边靠着顶点的情况
	repeat(i,0,a.size())a[i].x=l-a[i].x; // 水平翻折
	sort(a.begin(),a.end(),[](vec a,vec b){return a.x<b.x;});
	repeat(i,0,a.size())
		work(i);
	// 小矩形左右边都不靠顶点的情况
	sort(a.begin(),a.end(),[](vec a,vec b){return a.y<b.y;});
	repeat(i,0,(int)a.size()-1)
		ans=max(ans,(a[i+1].y-a[i].y)*l);
	return ans;
}

平面最近点对 using 分治

  • \(O(n\log n)\)
lf ans;
bool cmp_y(vec a,vec b){return a.y<b.y;}
void work(int l,int r){
	#define upd(x,y) {ans=min(ans,(x-y).len());}
	if(r-l<=4){
		repeat(i,l,r)
		repeat(j,i+1,r)
			upd(a[i],a[j]);
		sort(a+l,a+r,cmp_y);
		return;
	}
	int m=(l+r)/2;
	lf midx=a[m].x;
	work(l,m); work(m,r);
	static vec b[N];
	inplace_merge(a+l,a+m,a+r,cmp_y);
	int t=0;
	repeat(i,l,r)
	if(abs(a[i].x-midx)<ans){
		repeat_back(j,0,t){
			if(a[i].y-b[j].y>ans)break;
			upd(a[i],b[j]);
		}
		b[t++]=a[i];
	}
}
lf nearest(int n){
	ans=1e20;
	sort(a,a+n,[](vec a,vec b){return a.x<b.x;});
	work(0,n);
	return ans;
}

最小圆覆盖 using 随机增量法

  • eps可能要非常小。随机化，均摊 \(O(n)\)
struct cir{ // 圆（结构体）
	vec v; lf r;
	bool out(vec b){ // 点a在圆外
		return (v-b).len()>r+eps;
	}
	cir(vec a){v=a; r=0;}
	cir(vec a,vec b){v=(a+b)*0.5; r=(v-a).len();}
	cir(vec a,vec b,vec c){ // 三个点的外接圆
		b=b-a,c=c-a;
		vec s=vec(b.sqr(),c.sqr())*0.5;
		lf d=1/cross(b,c);
		v=a+vec(s.x*c.y-s.y*b.y,s.y*b.x-s.x*c.x)*d;
		r=(v-a).len();
	}
};
cir RIA(vec a[],int n){
	repeat_back(i,2,n)swap(a[rand()%i],a[i]); // random_shuffle(a,a+n);
	cir c=cir(a[0]);
	repeat(i,1,n)if(c.out(a[i])){
		c=cir(a[i]);
		repeat(j,0,i)if(c.out(a[j])){
			c=cir(a[i],a[j]);
			repeat(k,0,j)if(c.out(a[k]))
				c=cir(a[i],a[j],a[k]);
		}
	}
	return c;
}

半面交 using S&I 算法

  • 编号从 0 开始，\(O(n\log n)\)
vector<vec> ans; // ans: output, shows a convex hull
namespace half{
line a[N]; int n; // (a[],n): input, the final area will be the left of the lines
deque<line> q;
void solve(){
	a[n++]=line(vec(inf,inf),vec(-inf,inf));
	a[n++]=line(vec(-inf,inf),vec(-inf,-inf));
	a[n++]=line(vec(-inf,-inf),vec(inf,-inf));
	a[n++]=line(vec(inf,-inf),vec(inf,inf));
	sort(a,a+n,[](line a,line b){
		if(a.th<b.th-eps)return 1;
		if(a.th<b.th+eps && b.contain(a.p1)==1)return 1;
		return 0;
	});
	n=unique(a,a+n,[](line a,line b){return abs(a.th-b.th)<eps;})-a;
	q.clear();
	#define r q.rbegin()
	repeat(i,0,n){
		while(q.size()>1 && !a[i].contain(r[0].PI(r[1])))q.pop_back();
		while(q.size()>1 && !a[i].contain(q[0].PI(q[1])))q.pop_front();
		q.push_back(a[i]);
	}
	while(q.size()>1 && !q[0].contain(r[0].PI(r[1])))q.pop_back();
	while(q.size()>1 && !r[0].contain(q[0].PI(q[1])))q.pop_front();
	#undef r
	ans.clear();
	repeat(i,0,(int)q.size()-1)ans<<q[i].PI(q[i+1]);
	ans<<q[0].PI(q.back());
}
}

struct of 整点直线

struct line{
	ll up,down,dx,dy; // y=(dy/dx)x+(up/down) or x=(up/down)
	void adjust(ll &x,ll &y){
		if(x<0)x=-x,y=-y;
		if(x==0)y=1;
		else if(y==0)x=1;
		else{
			ll d=abs(__gcd(x,y));
			x/=d; y/=d;
		}
	}
	line(ll x1,ll y1,ll x2,ll y2){
		dx=(x1-x2),dy=(y1-y2);
		adjust(dx,dy);
		if(dx!=0){
			up=-dy*x1+dx*y1;
			down=dx;
			adjust(up,down);
		}
		else{
			up=-dx*y1+dy*x1;
			down=dy;
			adjust(up,down);
		}
	}
	pii d(){return {dx,dy};} // 斜率
	pii d2(){ // 垂线斜率
		ll ddx=-dy,ddy=dx;
		adjust(ddx,ddy);
		return {ddx,ddy};
	}
	bool operator==(line b)const{
		return make_tuple(up,down,dx,dy)
			== make_tuple(b.up,b.down,b.dx,b.dy);
	}
};
struct h{ // Hash
	ll operator()(line a)const{
		return a.up+a.down*10000+a.dx*100000000+a.dy*1000000000000;
	}
};

整点向量线性基

  • n 个向量的集合 \(\{(x_i,y_i)\}\) 可以构造线性等价的两个向量的集合 \(\{(a_1,b_1),(a_2,b_2)\},(b_2=0)\)，即 \(\displaystyle\{\sum_{i=1}^n t_i(x_i,y_i)\mid t\in \Z^n\}=\{t_1(a_1,b_1)+t_2(a_2,b_2)\mid t\in \Z^2\}\)
  • linear::push(x,y): 添加向量
  • linear::query(x,y): 询问向量能否被线性表示
struct linear{
	ll a1,b1,a2; // (a1,b1),(a2,0)
	linear(){a1=b1=a2=0;}
	void push(ll x,ll y){
		ll A,B,d;
		exgcd(y,b1,d,A,B);
		a2=__gcd(a2,abs(y*a1-x*b1)/d);
		b1=d;
		a1=A*x+B*a1;
		if(a2)a1=(a1%a2+a2)%a2;
	}
	bool query(ll x,ll y){
		if(b1!=0 && y%b1==0){
			ll r=x-y/b1*a1;
			return (a2!=0 && r%a2==0) || a2==r;
		}
		else if(b1==0)
			return (a1!=0 && x%a1==0) || a1==x;
		else return false;
	}
};

曼哈顿最小生成树

  • input: n,a[i].x,a[i].y，a[i].p 是没用的。编号从 1 开始，\(O(n\log n)\)
DSU d;
int n,w[N],c[N];
struct node{
	int x,y,p;
}a[N],b[N];
vector<node> e;
int dist(int x,int y){
	return abs(a[x].x-a[y].x)+abs(a[x].y-a[y].y);
}
#define lb(x) (x&-x)
struct BIT{ // special
	int t[N];
	void init(){
		fill(t,t+n+1,0);
	}
	void insert(int x,int p){
		for(;x<=n;x+=lb(x))
		if(w[p]<=w[t[x]])
			t[x]=p;
	}
	int query(int x){
		int ans=0;
		for(;x!=0;x-=lb(x))
		if(w[t[x]]<=w[ans])
			ans=t[x];
		return ans;
	}
}bit; 
void work(){
	bit.init();
	repeat(i,1,n+1)c[i]=b[i].y; sort(c+1,c+n+1);
	sort(b+1,b+n+1,[](node a,node b){
		return pii(a.x,a.y)<pii(b.x,b.y);
	});
	repeat(i,1,n+1){
		int u=upper_bound(c+1,c+n+1,b[i].y)-c,j=bit.query(u);
		if(j)e.push_back({b[i].p,j,dist(b[i].p,j)});
		bit.insert(u,b[i].p);
	}
}
ll mmst(){
	w[0]=inf; e.clear(); d.init(n);
	repeat(i,1,n+1){
		b[i]={-a[i].x,a[i].x-a[i].y,i};
		w[i]=a[i].x+a[i].y;
	}
	work();
	repeat(i,1,n+1){
		b[i]={-a[i].y,a[i].y-a[i].x,i};
	}
	work();
	repeat(i,1,n+1){
		b[i]={a[i].y,-a[i].x-a[i].y,i};
		w[i]=a[i].x-a[i].y;
	}
	work();
	repeat(i,1,n+1){
		b[i]={-a[i].x,a[i].y+a[i].x,i};
	}
	work();
	sort(e.begin(),e.end(),[](node a,node b){
		return a.p<b.p;
	});
	ll ans=0;
	for(auto i:e)
	if(d[i.x]!=d[i.y]){
		d[i.x]=d[i.y],ans+=i.p;
	}
	return ans;
}

圆的离散化

  • refer to CTSC 07 高逸涵
  • 若干圆，任意两圆不相切，求未被圆覆盖的闭合图形个数
  • 将圆的上下顶点和两两圆的交点的y作为事件，取相邻事件中点 \(e[i]\)，分析其状态，对相邻的 \(e[i]\) 用并查集判连通
// poj 1688 但是 wa 了
DSU d;
vector<lf> ovo,e;
vector<pair<lf,lf> > rec[N];
vector<int> lab[N]; int labcnt,ans;
void segunion(vector<pair<lf,lf> > &a){ // 区间合并至最简
	if(a.empty())return;
	sort(a.begin(),a.end()); int pre=0;
	repeat(i,0,a.size()){
		if(a[i].fi>a[pre].se-eps)a[++pre]=a[i];
		else a[pre].se=max(a[pre].se,a[i].se);
	}
	a.erase(a.begin()+pre+1,a.end());
}
void segcomplement(vector<pair<lf,lf> > &a){ // 区间取反
	a.push_back(pair<lf,lf>(0,inf));
	repeat_back(i,0,a.size()-1){
		a[i+1].fi=a[i].se;
		a[i].se=a[i].fi;
	}
	a[0].fi=-inf;
}
void Solve(){
	int n=read(); ovo.clear(); e.clear(); labcnt=ans=0;
	repeat(i,0,n){
		a[i].v.x=read(),a[i].v.y=read(),a[i].r=read();
		ovo.push_back(a[i].v.y-a[i].r);
		ovo.push_back(a[i].v.y+a[i].r);
	}
	repeat(i,0,n)
	repeat(j,i+1,n)
	if((a[i].v-a[j].v).len()<a[i].r+a[j].r){
		vec A,B; a[i].PI(a[j],A,B);
		ovo.push_back(A.y);
		ovo.push_back(B.y);
	}
	sort(ovo.begin(),ovo.end());
	e.push_back(-inf);
	repeat(i,0,ovo.size()-1)
		e.push_back((ovo[i]+ovo[i+1])/2);
	e.push_back(inf);
	repeat(j,0,e.size()){
		rec[j].clear();
		repeat(i,0,n)
		if(abs(a[i].v.y-e[j])<a[i].r-eps){
			lf d=sqrt(a[i].r*a[i].r-(a[i].v.y-e[j])*(a[i].v.y-e[j]));
			rec[j].push_back(pair<lf,lf>(a[i].v.x-d,a[i].v.x+d));
		}
		segunion(rec[j]); 
		segcomplement(rec[j]);
		lab[j].assign(rec[j].size(),0);
		repeat(i,0,lab[j].size())lab[j][i]=labcnt++;
	}
	d.init(labcnt);
	repeat(i,0,e.size()-1){
		unsigned p1=0,p2=0;
		while(p1<rec[i].size() && p2<rec[i+1].size()){
			if(rec[i][p1].se>rec[i+1][p2].fi-eps && rec[i][p1].fi<rec[i+1][p2].se+eps){
				int x=lab[i][p1],y=lab[i+1][p2];
				if(d[x]!=d[y]){
					ans++;
					d[x]=d[y];
				}
			}
			(rec[i][p1].se<rec[i+1][p2].se?p1:p2)++;
		}
	}
	printf("%d\n",labcnt-ans-1);
}

Delaunay 三角剖分

  • 编号从 0 开始，\(O(n\log n)\)
const lf eps=1e-8;
struct vec{
	lf x,y; int id;
	explicit vec(lf a=0,lf b=0,int c=-1):x(a),y(b),id(c){}
	bool operator<(const vec &a)const{
		return x<a.x || (abs(x-a.x)<eps && y<a.y);
	}
	bool operator==(const vec &a)const{
		return abs(x-a.x)<eps && abs(y-a.y)<eps;
	}
	lf dist2(const vec &b){
		return (x-b.x)*(x-b.x)+(y-b.y)*(y-b.y);
	}
};
struct vec3D{
	lf x,y,z;
	explicit vec3D(lf a=0,lf b=0,lf c=0):x(a),y(b),z(c){}
	vec3D(const vec &v){x=v.x,y=v.y,z=v.x*v.x+v.y*v.y;}
	vec3D operator-(const vec3D &a)const{
		return vec3D(x-a.x,y-a.y,z-a.z);
	}
};
struct edge{
	int id;
	list<edge>::iterator c;
	edge(int id=0){this->id=id;}
};
int cmp(lf v){return abs(v)>eps?(v>0?1:-1):0;}
lf cross(const vec &o,const vec &a,const vec &b){
	return(a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
lf dot(const vec3D &a,const vec3D &b){return a.x*b.x+a.y*b.y+a.z*b.z;}
vec3D cross(const vec3D &a,const vec3D &b){
	return vec3D(a.y*b.z-a.z*b.y,-a.x*b.z+a.z*b.x,a.x*b.y-a.y*b.x);
}
vector<pii> ans; // 三角剖分结果
struct DT{ // 使用方法：直接solve()
	list<edge> a[N]; vec v[N]; int n;
	void solve(int _n,vec _v[]){
		n=_n;
		copy(_v,_v+n,v);
		sort(v,v+n);
		divide(0,n-1);
		ans.clear();
		for(int i=0;i<n;i++){
			for(auto p:a[i]){
				if(p.id<i)continue;
				ans.push_back({v[i].id,v[p.id].id});
			}
		}
	}
	int incircle(const vec &a,vec b,vec c,const vec &v){
		if(cross(a,b,c)<0)swap(b,c);
		vec3D a3(a),b3(b),c3(c),p3(v);
		b3=b3-a3,c3=c3-a3,p3=p3-a3;
		vec3D f=cross(b3,c3);
		return cmp(dot(p3,f));
	}
	int intersection(const vec &a,const vec &b,const vec &c,const vec &d){
		return cmp(cross(a,c,b))*cmp(cross(a,b,d))>0 &&
			cmp(cross(c,a,d))*cmp(cross(c,d,b))>0;
	}
	void addedge(int u,int v){
		a[u].push_front(edge(v));
		a[v].push_front(edge(u));
		a[u].begin()->c=a[v].begin();
		a[v].begin()->c=a[u].begin();
	}
	void divide(int l,int r){
		if(r-l<=2){
			for(int i=l;i<=r;i++)
				for(int j=i+1;j<=r;j++)addedge(i,j);
			return;
		}
		int mid=(l+r)/2;
		divide(l,mid); divide(mid+1,r);
		int nowl=l,nowr=r;
		for(int update=1;update;){
			update=0;
			vec vl=v[nowl],vr=v[nowr];
			for(auto i:a[nowl]){
				vec t=v[i.id];
				lf v=cross(vr,vl,t);
				if(cmp(v)>0 || (cmp(v)== 0 && vr.dist2(t)<vr.dist2(vl))){
					nowl=i.id,update=1;
					break;
				}
			}
			if(update)continue;
			for(auto i:a[nowr]){
				vec t=v[i.id];
				lf v=cross(vl,vr,t);
				if(cmp(v)<0 || (cmp(v)== 0 && vl.dist2(t)<vl.dist2(vr))){
					nowr=i.id,update=1;
					break;
				}
			}
		}
		addedge(nowl,nowr);
		while(1){
			vec vl=v[nowl],vr=v[nowr];
			int ch=-1,side=0;
			for(auto i:a[nowl])
			if(cmp(cross(vl,vr,v[i.id]))>0 && (ch==-1 || incircle(vl,vr,v[ch],v[i.id])<0)){
				ch=i.id,side=-1;
			}
			for(auto i:a[nowr])
			if(cmp(cross(vr,v[i.id],vl))>0 && (ch==-1 || incircle(vl,vr,v[ch],v[i.id])<0)){
				ch=i.id,side=1;
			}
			if(ch==-1)break;
			if(side==-1){
				for(auto it=a[nowl].begin();it!=a[nowl].end();){
					if(intersection(vl,v[it->id],vr,v[ch])){
						a[it->id].erase(it->c);
						a[nowl].erase(it++);
					}
					else it++;
				}
				nowl=ch;
				addedge(nowl,nowr);
			}
			else{
				for(auto it=a[nowr].begin();it!=a[nowr].end();){
					if(intersection(vr,v[it->id],vl,v[ch])){
						a[it->id].erase(it->c);
						a[nowr].erase(it++);
					}
					else it++;
				}
				nowr=ch;
				addedge(nowl,nowr);
			}
		}
	}
}dt;
  • 可以求最小生成树
vec a[N]; DSU d;
vector<int> e[N]; // 最小生成树结果
void mst(){ // 求最小生成树
	dt.solve(n,a);
	sort(ans.begin(),ans.end(),[](const pii &A,const pii &B){
		return a[A.fi].dist2(a[A.se])<a[B.fi].dist2(a[B.se]);
	});
	d.init(n);
	for(auto i:ans)
	if(d[i.fi]!=d[i.se]){
		e[i.fi].push_back(i.se);
		e[i.se].push_back(i.fi);
		d[i.fi]=d[i.se];
	}
}

struct of 三维向量

  • trunc(K) 返回 K 在 *this 上的投影向量
  • rotate(P,L,th) 返回点 P 绕轴 (O,L) 旋转 th 弧度后的点
  • rotate(P,L0,L1,th) 返回点 P 绕轴 (L0,L1) 旋转 th 弧度后的点
struct vec{
	lf x,y,z; vec(){} vec(lf x,lf y,lf z):x(x),y(y),z(z){}
	vec operator-(vec b){return vec(x-b.x,y-b.y,z-b.z);}
	vec operator+(vec b){return vec(x+b.x,y+b.y,z+b.z);}
	vec operator*(lf k){return vec(k*x,k*y,k*z);}
	bool operator<(vec b)const{return make_tuple(x,y,z)<make_tuple(b.x,b.y,b.z);}
	lf sqr(){return x*x+y*y+z*z;}
	lf len(){return sqrt(x*x+y*y+z*z);}
	vec trunc(lf k=1){return *this*(k/len());}
	vec trunc(vec k){return *this*(dot(*this,k)/sqr());}
	friend vec cross(vec a,vec b){
		return vec(
			a.y*b.z-a.z*b.y,
			a.z*b.x-a.x*b.z,
			a.x*b.y-a.y*b.x
		);
	}
	friend lf dot(vec a,vec b){return a.x*b.x+a.y*b.y+a.z*b.z;}
	friend vec rotate(vec p,vec l,lf th){
		struct four{
			lf r; vec v;
			four operator*(four b){
				return {r*b.r-dot(v,b.v),v*b.r+b.v*r+cross(v,b.v)};
			}
		};
		l=l.trunc();
		four P={0,p};
		four Q1={cos(th/2),l*sin(th/2)};
		four Q2={cos(th/2),vec()-l*sin(th/2)};
		return ((Q1*P)*Q2).v;
	}
	friend vec rotate(vec p,vec l0,vec l1,lf th){
		return rotate(p-l0,l1-l0,th)+l0;
	}
	void print(){printf("%.12f %.12f %.12f\n",x,y,z);}
};

三维凸包

  • 将所有凸包上的面放入面集 f 中，其中 face::p[i] 作为 a 的下标，\(O(n^2)\)
const lf eps=1e-9;
struct vec{
	lf x,y,z;
	vec(lf x=0,lf y=0,lf z=0):x(x),y(y),z(z){};
	vec operator-(vec b){return vec(x-b.x,y-b.y,z-b.z);}
	lf len(){return sqrt(x*x+y*y+z*z);}
	void shake(){ // 微小扰动
		x+=(rand()*1.0/RAND_MAX-0.5)*eps;
		y+=(rand()*1.0/RAND_MAX-0.5)*eps;
		z+=(rand()*1.0/RAND_MAX-0.5)*eps;
	}
}a[N];
vec cross(vec a,vec b){
	return vec(
		a.y*b.z-a.z*b.y,
		a.z*b.x-a.x*b.z,
		a.x*b.y-a.y*b.x);
}
lf dot(vec a,vec b){return a.x*b.x+a.y*b.y+a.z*b.z;}
struct face{
	int p[3];
	vec normal(){ // 法向量
		return cross(a[p[1]]-a[p[0]],a[p[2]]-a[p[0]]);
	}
	lf area(){return normal().len()/2.0;}
};
vector<face> f;
bool see(face f,vec v){
	return dot(v-a[f.p[0]],f.normal())>0;
}
void convex(vec a[],int n){
	static vector<face> c;
	static bool vis[N][N];
	repeat(i,0,n)a[i].shake(); // 防止四点共面
	f.clear();
	f.push_back((face){0,1,2});
	f.push_back((face){0,2,1});
	repeat(i,3,n){
		c.clear();
		repeat(j,0,f.size()){
			bool t=see(f[j],a[i]);
			if(!t) // 加入背面
				c.push_back(f[j]);
			repeat(k,0,3){
				int x=f[j].p[k],y=f[j].p[(k+1)%3];
				vis[x][y]=t;
			}
		}
		repeat(j,0,f.size())
		repeat(k,0,3){
			int x=f[j].p[k],y=f[j].p[(k+1)%3];
			if(vis[x][y] && !vis[y][x]) // 加入新面
				c.push_back((face){x,y,i});
		}
		f.swap(c);
	}
}
posted @ 2021-07-07 01:16  axiomofchoice  阅读(211)  评论(0)  编辑  收藏  举报