bzoj 3053 HDU 4347 : The Closest M Points kd树

时间:2023-06-03 14:53:56

bzoj 3053 HDU 4347 : The Closest M Points  kd树

题目大意:求k维空间内某点的前k近的点。

就是一般的kd树,根据实测发现,kd树的两种建树方式,即按照方差较大的维度分开(建树常数大)或者每一位轮换分割(询问常数大),后者更快也更好些,以后就果断写第二种了。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
using namespace std;
#define MAXN 510000
#define MAXT MAXN
#define MAXM 6
#define sqr(x) ((qword)(x)*(x))
#define INF 0x3f3f3f3f
typedef long long qword;
int n,m;
struct point
{
int a[MAXM];
qword dis;
void pm()
{
printf("(%d",a[]);
for (int i=;i<m;i++)
printf(",%d",a[i]);
printf(")");
}
void pm2()
{
for (int i=;i<m-;i++)
printf("%d ",a[i]);
printf("%d\n",a[m-]);
}
}pl[MAXN];
bool cmp_0(const point &p1,const point &p2){
return p1.a[]<p2.a[];
}
bool cmp_1(const point &p1,const point &p2){
return p1.a[]<p2.a[];
}
bool cmp_2(const point &p1,const point &p2){
return p1.a[]<p2.a[];
}
bool cmp_3(const point &p1,const point &p2){
return p1.a[]<p2.a[];
}
bool cmp_4(const point &p1,const point &p2){
return p1.a[]<p2.a[];
}
bool cmp_d(const point &p1,const point &p2){
return p1.dis<p2.dis;
}
struct kdt_node
{
int a[MAXM];
int dd,mxp[MAXM],mnp[MAXM];
int ch[];
int ptr;
}kdt[MAXT];
qword get_dis(point &pt,kdt_node &pn)
{
qword ret=;
for (int j=;j<m;j++)
if (pt.a[j]<pn.mnp[j] || pt.a[j]>pn.mxp[j])
ret+=min(sqr(pt.a[j]-pn.mxp[j]),sqr(pt.a[j]-pn.mnp[j]));
return ret;
}
qword get_dis(point &p1,point &p2)
{
qword ret=;
for (int i=;i<m;i++)
ret+=sqr(p1.a[i]-p2.a[i]);
return ret;
}
int topt=;
void Build_kdt(int &now,int l,int r,int d)
{
if (l>r)return;
now=++topt;
int i,j;
for (j=;j<m;j++)kdt[now].mxp[j]=-INF,kdt[now].mnp[j]=INF;
for (i=l;i<=r;i++)
for (j=;j<m;j++)
{
kdt[now].mxp[j]=max(kdt[now].mxp[j],pl[i].a[j]);
kdt[now].mnp[j]=min(kdt[now].mnp[j],pl[i].a[j]);
}
/* double ave[MAXM];
double sqv[MAXM];
memset(ave,0,sizeof(ave));
memset(sqv,0,sizeof(sqv));
for (i=l;i<=r;i++)
for (j=0;j<m;j++)
ave[j]+=pl[i].a[j];
for (j=0;j<m;j++)
ave[j]/=(r-l+1);
for (i=l;i<=r;i++)
for (j=0;j<m;j++)
sqv[j]+=sqr(pl[i].a[j]-ave[j]);
kdt[now].dd=0;
for (j=0;j<m;j++)
if (sqv[j]>sqv[kdt[now].dd])
kdt[now].dd=j;*/
kdt[now].dd=d;
switch (kdt[now].dd)
{
case :nth_element(&pl[l],&pl[(l+r)>>],&pl[r+],cmp_0);break;
case :nth_element(&pl[l],&pl[(l+r)>>],&pl[r+],cmp_1);break;
case :nth_element(&pl[l],&pl[(l+r)>>],&pl[r+],cmp_2);break;
case :nth_element(&pl[l],&pl[(l+r)>>],&pl[r+],cmp_3);break;
case :nth_element(&pl[l],&pl[(l+r)>>],&pl[r+],cmp_4);break;
}
kdt[now].ptr=(l+r)>>;
Build_kdt(kdt[now].ch[],l,((r+l)>>)-,(d+)%m);
Build_kdt(kdt[now].ch[],((r+l)>>)+,r,(d+)%m);
}
point h[MAXN];
int atot;
int toph=;
void search_point(int now,point &pt)
{
if (!now)return ;
qword cdis=get_dis(pt,pl[kdt[now].ptr]);
if (toph<atot || cdis<h[].dis)
{
if (toph==atot)pop_heap(&h[],&h[toph--],cmp_d);
h[toph]=pl[kdt[now].ptr];
h[toph].dis=cdis;
push_heap(&h[],&h[++toph],cmp_d);
}
int t;
if (get_dis(pt,kdt[kdt[now].ch[]]) < get_dis(pt,kdt[kdt[now].ch[]]))
t=;
else
t=;
search_point(kdt[now].ch[t],pt);
if (toph<atot || get_dis(pt,kdt[kdt[now].ch[t^]]) < h[].dis)
{
search_point(kdt[now].ch[t^],pt);
}
}
int main()
{
// freopen("input.txt","r",stdin);
//freopen("output.txt","w",stdout);
int i,j,k,x,y,z;
int root=;
while (~scanf("%d%d",&n,&m))
{
for (i=;i<n;i++)
for (j=;j<m;j++)
scanf("%d",&pl[i].a[j]);
Build_kdt(root,,n-,);
int q;
scanf("%d",&q);
point pt;
for (i=;i<q;i++)
{
for (j=;j<m;j++)
scanf("%d",&pt.a[j]);
scanf("%d",&atot);
search_point(root,pt);
printf("the closest %d points are:\n",atot);
while (toph)
{
pop_heap(&h[],&h[toph--],cmp_d);
}
for (j=;j<atot;j++)
h[j].pm2();
}
}
}