当前位置：首页 > article >正文

WSPD：平面最近邻+t-spanner+近似欧氏距离MST（程设实习）

article 2024/12/29 10:09:46

~~哎呀打字麻烦死了还是贴PPT吧~~

前言

感觉这个东西还是有点厉害的。

定义

在这里插入图片描述

构造

在这里插入图片描述
可以证明这样构造的WSPD大小为 $O(n\epsilon^{-d}\log\triangle)$ ， $\triangle$ 为值域。

平面最近邻

构造 $\epsilon=1$ 的WASD，只需对于所有 $A_i|=|B_i|=1$ 的里面的两个点的距离向答案贡献一次即可。

代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define uint unsigned int
#define debug(...) fprintf(stderr,__VA_ARGS__)
#define ok debug("ok\n")
inline ll read(){
	ll x(0),f(1);char c=getchar();
	while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
	while(isdigit(c)){x=(x<<1)+(x<<3)+c-'0';c=getchar();}
	return f*x;
}

const int N=2e5+100;
const double eps=1.0;

double mx=1<<30;
int n,m;

struct point{
	double x[3];
}; 
void print(point o,int flag){
	printf("(%.0lf %.0lf %.0lf) ",o.x[0],o.x[1],o.x[2]);
	if(flag) puts("");
}
inline double dis(const point &u,const point &v){
	double res=0;
	for(int i=0;i<3;i++){
		res+=(u.x[i]-v.x[i])*(u.x[i]-v.x[i]);
	}
	return sqrt(res);
}
struct tree{
	double l[3],r[3];
	point o;
	int son[8];
}tr[N*20];
const double c=sqrt(3);
inline double diam(int x){
	if(tr[x].son[0]==-1) return 0;
	return c*(tr[x].r[0]-tr[x].l[0]);
}
inline double dis(const tree &u,const tree &v){
	
	double ans=0;
	for(int i=0;i<3;i++){
		double d;
		if(u.r[i]>v.l[i]&&v.r[i]>u.l[i]) d=0;
		else if(u.l[i]>v.r[i]) d=u.l[i]-v.r[i];
		else d=v.l[i]-u.r[i];
		ans+=d*d;
		//printf("i=%d d=%lf flag=%d %d %lf>%lf\n",i,d,u.r[i]+1e-9>v.l[i],u.r[i]<v.r[i]+1e-9u.r[i],v.r[i]+1e-9);
	}
	return sqrt(ans);
}

int tot;
int build(vector<point>ve,vector<double> le,vector<double> ri){
	if(ve.empty()) return 0;
	int now=++tot;
	for(int i=0;i<3;i++) tr[now].l[i]=le[i],tr[now].r[i]=ri[i];
	if(ve.size()==1){
		tr[now].o=ve[0];
		tr[now].son[0]=-1;
		return now;
	}
	vector<double> mid;
	for(int i=0;i<3;i++) mid.push_back((le[i]+ri[i])/2);
	vector<point>v[8];
	for(point o:ve){
		int id=0;
		for(int i=0;i<3;i++){
			if(o.x[i]>mid[i]) id|=(1<<i);
		}
		v[id].push_back(o);
	}
	for(int id=0;id<8;id++){
		vector<double>l,r;
		for(int i=0;i<3;i++){
			if(id>>i&1){
				l.push_back(mid[i]);
				r.push_back(ri[i]);
			}
			else{
				l.push_back(le[i]);
				r.push_back(mid[i]);
			}
		}
		tr[now].son[id]=build(v[id],l,r);
	}
	return now;
}


double ans=1e18;
vector<point>tmp;
void dfs(int x){
	if(!x) return;
	if(tr[x].son[0]==-1){
		tmp.push_back(tr[x].o);
		return;
	}
	for(int i=0;i<8;i++) dfs(tr[x].son[i]);
}
void getpt(int x){
	tmp.clear();
	dfs(x);
	for(point o:tmp) print(o,0);
	puts("");
}
inline void get(int x,int y){
	//printf("\nget: (%d %d)\n",x,y);
	//getpt(x);
	//getpt(y);
	//puts("");
	if(tr[x].son[0]==-1&&tr[y].son[0]==-1){
		ans=min(ans,dis(tr[x].o,tr[y].o));
		//print(tr[x].o,0);
		//print(tr[y].o,1);
	}
	return;
}
void find(int u,int v){
	if(u*v==0) return;
	if(u==v&&tr[u].son[0]==-1) return;
	if(diam(u)<diam(v)) swap(u,v);
	//printf("find: (%d %d) diam=%lf %lf dis=%lf\n",u,v,diam(u),diam(v),dis(tr[u],tr[v]));
	//getpt(u);
	//getpt(v);
	if(diam(u)<=eps*dis(tr[u],tr[v])){//find a WSPD
		get(u,v);
		return;
	}
	for(int i=0;i<8;i++) find(tr[u].son[i],v);
}

int main() {
	
	//freopen("a.in","r",stdin);
	//freopen("a.out","w",stdout);
    	
	n=read();
	vector<point>ve;
	ve.resize(n);
	for(point &o:ve){
		for(int i=0;i<3;i++) o.x[i]=read();
	}
	
	if(0)for(int i=0;i<n;i++){
		for(int j=i+1;j<n;j++) printf("(%d %d) dis=%lf\n",i,j,dis(ve[i],ve[j]));
	}
	
	vector<double>l,r;
	for(int i=0;i<3;i++){
		l.push_back(0);
		r.push_back(mx);
	}
	int rt=build(ve,l,r);
	find(rt,rt);
    
    printf("%.2lf\n",ans);
    return 0;
}