三维空间曼哈顿最小生成树

题面
二维空间最小生成树的解决方法是以每个点为中心将二维平面分成米字型的八个区域，在每个区域里求出距离最近的点连边。这样做的原理是每个区域满足中心点与更远的点连边不优。
同样地，将三维空间每个卦限划分为四个三棱锥型区域，中心点在每个区域中找最近点连边。划分方式见下图：

图中1,2,3,4为四个三棱锥区域。
每个三棱锥区域都可以用三个面来描述，因此可以转为三维偏序限制。对于三维偏序问题，用CDQ分治即可。
时间复杂度 $O(n{\log }_2^{2}n)$，空间复杂度 $O(n)$

#include<stdio.h>
#include<vector>
#include<algorithm>
using namespace std;
#define R register int
#define I inline
#define N 50001
#define INF -999999999
int X[N],Y[N],Z[N],Tx[N],Ty[N],Tz[N],f[N];
struct Edge{
	int Start,End,Len;
	I friend bool operator<(Edge A,Edge B){
		return A.Len<B.Len;
	}
};
vector<Edge>E;
I void AddEdge(int x,int y){
	Edge res;
	res.Start=x;
	res.End=y;
	int a=X[x]-X[y],b=Y[x]-Y[y],c=Z[x]-Z[y];
	res.Len=(a<0?-a:a)+(b<0?-b:b)+(c<0?-c:c);
	E.push_back(res);
}
I bool Compare(int x,int y){
	if(Tx[x]!=Tx[y]){
		return Tx[x]<Tx[y];
	}
	if(Ty[x]!=Ty[y]){
		return Ty[x]<Ty[y];
	}
	return Tz[x]<Tz[y];
}
class CDQ{
	int dist[N],c[N],d[N],id[N],tem[N],cur[N],curD[N],Tot;
	I void Insert(const int q){
		int v=X[q]+Y[q]+Z[q];
		for(R i=Tz[q];i<=Tot;i+=i&-i){
			if(v>c[i]){
				c[i]=v;
				d[i]=q;
			}
		}
	}
	I int GetAns(int x){
		int res=INF,resid=-1;
		for(R i=x;i!=0;i&=i-1){
			if(c[i]>res){
				res=c[i];
				resid=d[i];
			}
		}
		return resid;
	}
	I void Clear(int x){
		for(R i=x;i<=Tot;i+=i&-i){
			c[i]=INF;
			d[i]=-1;
		}
	}
	I void MergeSort(int l,int r){
		if(l!=r){
			int mid=l+r>>1,p1=l,p2;
			MergeSort(l,mid);
			p2=mid+1;
			MergeSort(p2,r);
			for(R i=l;i<=r;i++){
				if(p1>mid||p2<=r&&Ty[id[p1]]>Ty[id[p2]]){
					int x=GetAns(Tz[id[p2]]);
					if(x!=-1){
						int y=X[x]+Y[x]+Z[x];
						if(y>curD[id[p2]]){
							curD[id[p2]]=y;
							cur[id[p2]]=x;
						}
					}
					tem[i]=id[p2];
					p2++;
				}else{
					Insert(id[p1]);
					tem[i]=id[p1];
					p1++;
				}
			}
			for(R i=l;i<=mid;i++){
				Clear(Tz[id[i]]);
			}
			for(R i=l;i<=r;i++){
				id[i]=tem[i];
			}
		}
	}
	public:
		I void Solve(const int n){
			for(R i=0;i!=n;i++){
				dist[i]=Tz[i];
				id[i]=i;
				curD[i]=INF;
				cur[i]=-1;
			}
			sort(dist,dist+n);
			Tot=unique(dist,dist+n)-dist;
			for(R i=0;i!=n;i++){
				Tz[i]=upper_bound(dist,dist+Tot,Tz[i])-dist;
			}
			for(R i=1;i<=Tot;i++){
				c[i]=INF;
				d[i]=-1;
			}
			sort(id,id+n,Compare);
			MergeSort(0,n-1);
			for(R i=0;i!=n;i++){
				if(cur[i]!=-1){
					AddEdge(cur[i],i);
				}
			}
		}
}C;
I void RotateXY(const int n){
	for(R i=0;i!=n;i++){
		int t=Y[i];
		Y[i]=X[i];
		X[i]=-t;
	}
}
I void RotateXYZ(const int n){
	for(R i=0;i!=n;i++){
		int t=X[i];
		X[i]=Z[i];
		Z[i]=Y[i];
		Y[i]=t;
	}
}
I void Solve(const int n){
	for(R i=0;i!=n;i++){
		Tx[i]=Y[i]+Z[i]-X[i];
		Ty[i]=X[i]+Z[i]-Y[i];
		Tz[i]=X[i]+Y[i]-Z[i];
	}
	C.Solve(n);
	for(R i=0;i!=3;i++){
		for(R j=0;j!=n;j++){
			Tx[j]=X[j]-Y[j]-Z[j];
			Ty[j]=Y[j];
			Tz[j]=Z[j];
		}
		C.Solve(n);
		RotateXYZ(n);
	}
}
I int GetF(int x){
	if(f[x]==x){
		return x;
	}
	f[x]=GetF(f[x]);
	return f[x];
}
int main(){
	int n;
	scanf("%d",&n);
	for(R i=0;i!=n;i++){
		scanf("%d%d%d",X+i,Y+i,Z+i);
		f[i]=i;
	}
	Solve(n);
	for(R i=0;i!=3;i++){
		RotateXY(n);
		Solve(n);
	}
	sort(E.begin(),E.end());
	long long ans=0;
	for(vector<Edge>::iterator T=E.begin();T!=E.end();T++){
		if(GetF(T->Start)!=GetF(T->End)){
			ans+=T->Len;
			f[f[T->Start]]=T->End;
		}
	}
	printf("%lld",ans);
	return 0;
}