轉(zhuǎn)載:這個(gè)問(wèn)題很容易理解,似乎也不難解決。我們只要將每一點(diǎn)與其他n-1個(gè)點(diǎn)的距離算出,找出達(dá)到最小距離的兩個(gè)點(diǎn)即可。然而,這樣做效率太低,需要
O(n
2)的計(jì)算時(shí)間。在
問(wèn)題的計(jì)算復(fù)雜性中我們可以看到,該問(wèn)題的計(jì)算時(shí)間下界為
Ω(
nlog
n)。這個(gè)下界引導(dǎo)我們?nèi)フ覇?wèn)題的一個(gè)
θ(
nlog
n)算法。
這個(gè)問(wèn)題顯然滿(mǎn)足分治法的第一個(gè)和第二個(gè)適用條件,我們考慮將所給的平面上n個(gè)點(diǎn)的集合S分成2個(gè)子集S1和S2,每個(gè)子集中約有n/2個(gè)點(diǎn),·然后在每個(gè)子集中遞歸地求其最接近的點(diǎn)對(duì)。在這里,一個(gè)關(guān)鍵的問(wèn)題是如何實(shí)現(xiàn)分治法中的合并步驟,即由S1和S2的最接近點(diǎn)對(duì),如何求得原集合S中的最接近點(diǎn)對(duì),因?yàn)镾1和S2的最接近點(diǎn)對(duì)未必就是S的最接近點(diǎn)對(duì)。如果組成S的最接近點(diǎn)對(duì)的2個(gè)點(diǎn)都在S1中或都在S2中,則問(wèn)題很容易解決。但是,如果這2個(gè)點(diǎn)分別在S1和S2中,則對(duì)于S1中任一點(diǎn)p,S2中最多只有n/2個(gè)點(diǎn)與它構(gòu)成最接近點(diǎn)對(duì)的候選者,仍需做n2/4次計(jì)算和比較才能確定S的最接近點(diǎn)對(duì)。因此,依此思路,合并步驟耗時(shí)為O(n2)。整個(gè)算法所需計(jì)算時(shí)間T(n)應(yīng)滿(mǎn)足:
T(n)=2T(n/2)+O(n2)
它的解為T(mén)(n)=O(n2),即與合并步驟的耗時(shí)同階,顯示不出比用窮舉的方法好。從解遞歸方程的套用公式法,我們看到問(wèn)題出在合并步驟耗時(shí)太多。這啟發(fā)我們把注意力放在合并步驟上。
為了使問(wèn)題易于理解和分析,我們先來(lái)考慮一維的情形。此時(shí)S中的n個(gè)點(diǎn)退化為x軸上的n個(gè)實(shí)數(shù)x1,x2,..,xn。最接近點(diǎn)對(duì)即為這n個(gè)實(shí)數(shù)中相差最小的2個(gè)實(shí)數(shù)。我們顯然可以先將x1,x2,..,xn排好序,然后,用一次線(xiàn)性?huà)呙杈涂梢哉页鲎罱咏c(diǎn)對(duì)。這種方法主要計(jì)算時(shí)間花在排序上,因此如在排序算法中所證明的,耗時(shí)為O(nlogn)。然而這種方法無(wú)法直接推廣到二維的情形。因此,對(duì)這種一維的簡(jiǎn)單情形,我們還是嘗試用分治法來(lái)求解,并希望能推廣到二維的情形。
假設(shè)我們用x軸上某個(gè)點(diǎn)m將S劃分為2個(gè)子集S1和S2,使得S1={x∈S|x≤m};S2={x∈S|x>m}。這樣一來(lái),對(duì)于所有p∈S1和q∈S2有p<q。
遞歸地在S1和S2上找出其最接近點(diǎn)對(duì){p1,p2}和{q1,q2},并設(shè)δ=min{|p1-p2|,|q1-q2|},S中的最接近點(diǎn)對(duì)或者是{p1,p2},或者是{q1,q2},或者是某個(gè){p3,q3},其中p3∈S1且q3∈S2。如圖1所示。

圖1 一維情形的分治法
我們注意到,如果S的最接近點(diǎn)對(duì)是{p3,q3},即|p3-q3|<δ,則p3和q3兩者與m的距離不超過(guò)δ,即|p3-m|<δ,|q3-m|<δ,也就是說(shuō),p3∈(m-δ,m],q3∈(m,m+δ]。由于在S1中,每個(gè)長(zhǎng)度為δ的半閉區(qū)間至多包含一個(gè)點(diǎn)(否則必有兩點(diǎn)距離小于δ),并且m是S1和S2的分割點(diǎn),因此(m-δ,m]中至多包含S中的一個(gè)點(diǎn)。同理,(m,m+δ]中也至多包含S中的一個(gè)點(diǎn)。由圖1可以看出,如果(m-δ,m]中有S中的點(diǎn),則此點(diǎn)就是S1中最大點(diǎn)。同理,如果(m,m+δ]中有S中的點(diǎn),則此點(diǎn)就是S2中最小點(diǎn)。因此,我們用線(xiàn)性時(shí)間就能找到區(qū)間(m-δ,m]和(m,m+δ]中所有點(diǎn),即p3和q3。從而我們用線(xiàn)性時(shí)間就可以將S1的解和S2的解合并成為S的解。也就是說(shuō),按這種分治策略,合并步可在O(n)時(shí)間內(nèi)完成。這樣是否就可以得到一個(gè)有效的算法了呢?還有一個(gè)問(wèn)題需要認(rèn)真考慮,即分割點(diǎn)m的選取,及S1和S2的劃分。選取分割點(diǎn)m的一個(gè)基本要求是由此導(dǎo)出集合S的一個(gè)線(xiàn)性分割,即S=S1∪S2 ,S1∩S2=Φ,且S1
{x|x≤m};S2
{x|x>m}。容易看出,如果選取m=[max(S)+min(S)]/2,可以滿(mǎn)足線(xiàn)性分割的要求。選取分割點(diǎn)后,再用O(n)時(shí)間即可將S劃分成S1={x∈S|x≤m}和S2={x∈S|x>m}。然而,這樣選取分割點(diǎn)m,有可能造成劃分出的子集S1和S2的不平衡。例如在最壞情況下,|S1|=1,|S2|=n-1,由此產(chǎn)生的分治法在最壞情況下所需的計(jì)算時(shí)間T(n)應(yīng)滿(mǎn)足遞歸方程:
T(n)=T(n-1)+O(n)
它的解是T(n)=O(n2)。這種效率降低的現(xiàn)象可以通過(guò)分治法中“平衡子問(wèn)題”的方法加以解決。也就是說(shuō),我們可以通過(guò)適當(dāng)選擇分割點(diǎn)m,使S1和S2中有大致相等個(gè)數(shù)的點(diǎn)。自然地,我們會(huì)想到用S的n個(gè)點(diǎn)的坐標(biāo)的中位數(shù)來(lái)作分割點(diǎn)。在選擇算法中介紹的選取中位數(shù)的線(xiàn)性時(shí)間算法使我們可以在O(n)時(shí)間內(nèi)確定一個(gè)平衡的分割點(diǎn)m。
至此,我們可以設(shè)計(jì)出一個(gè)求一維點(diǎn)集S中最接近點(diǎn)對(duì)的距離的算法CPAIR1如下。
function CPAIR1(S);
begin
if |S|=2 then δ=|x[2]-x[1]| // x[1..n]存放的是S中n個(gè)點(diǎn)的坐標(biāo)
else if (|S|=1)
then δ:=∞
else begin
m:=S中各點(diǎn)的坐標(biāo)值的中位數(shù);
構(gòu)造S1和S2,使S1={x∈S|x≤m},S2={x∈S|x>m};
δ1:=CPAIRI(S1);
δ2:=CPAIRI(S2);
p:=max(S1);
q:=min(S2);
δ:=min(δ1,δ2,q-p);
end;
return(δ);
end;
由以上的分析可知,該算法的分割步驟和合并步驟總共耗時(shí)O(n)。因此,算法耗費(fèi)的計(jì)算時(shí)間T(n)滿(mǎn)足遞歸方程:

解此遞歸方程可得T(n)=O(nlogn)。
這個(gè)算法看上去比用排序加掃描的算法復(fù)雜,然而這個(gè)算法可以向二維推廣。
下面我們來(lái)考慮二維的情形。此時(shí)S中的點(diǎn)為平面上的點(diǎn),它們都有2個(gè)坐標(biāo)值x和y。為了將平面上點(diǎn)集S線(xiàn)性分割為大小大致相等的2個(gè)子集S1和S2,我們選取一垂直線(xiàn)l:x=m來(lái)作為分割直線(xiàn)。其中m為S中各點(diǎn)x坐標(biāo)的中位數(shù)。由此將S分割為S1={p∈S|px≤m}和S2={p∈S|px>m}。從而使S1和S2分別位于直線(xiàn)l的左側(cè)和右側(cè),且S=S1∪S2 。由于m是S中各點(diǎn)x坐標(biāo)值的中位數(shù),因此S1和S2中的點(diǎn)數(shù)大致相等。
遞歸地在S1和S2上解最接近點(diǎn)對(duì)問(wèn)題,我們分別得到S1和S2中的最小距離δ1和δ2。現(xiàn)設(shè)δ=min(δ1,δ1)。若S的最接近點(diǎn)對(duì)(p,q)之間的距離d(p,q)<δ則p和q必分屬于S1和S2。不妨設(shè)p∈S1,q∈S2。那么p和q距直線(xiàn)l的距離均小于δ。因此,我們?nèi)粲肞1和P2分別表示直線(xiàn)l的左邊和右邊的寬為δ的2個(gè)垂直長(zhǎng)條,則p∈S1,q∈S2,如圖2所示。

圖2 距直線(xiàn)l的距離小于δ的所有點(diǎn)
在一維的情形,距分割點(diǎn)距離為δ的2個(gè)區(qū)間(m-δ,m](m,m+δ]中最多各有S中一個(gè)點(diǎn)。因而這2點(diǎn)成為唯一的末檢查過(guò)的最接近點(diǎn)對(duì)候選者。二維的情形則要復(fù)雜些,此時(shí),P1中所有點(diǎn)與P2中所有點(diǎn)構(gòu)成的點(diǎn)對(duì)均為最接近點(diǎn)對(duì)的候選者。在最壞情況下有n2/4對(duì)這樣的候選者。但是P1和P2中的點(diǎn)具有以下的稀疏性質(zhì),它使我們不必檢查所有這n2/4對(duì)候選者。考慮P1中任意一點(diǎn)p,它若與P2中的點(diǎn)q構(gòu)成最接近點(diǎn)對(duì)的候選者,則必有d(p,q)<δ。滿(mǎn)足這個(gè)條件的P2中的點(diǎn)有多少個(gè)呢?容易看出這樣的點(diǎn)一定落在一個(gè)δ×2δ的矩形R中,如圖3所示。

圖3 包含點(diǎn)q的δ×2δ的矩形R
由δ的意義可知P2中任何2個(gè)S中的點(diǎn)的距離都不小于δ。由此可以推出矩形R中最多只有6個(gè)S中的點(diǎn)。事實(shí)上,我們可以將矩形R的長(zhǎng)為2δ的邊3等分,將它的長(zhǎng)為δ的邊2等分,由此導(dǎo)出6個(gè)(δ/2)×(2δ/3)的矩形。如圖4(a)所示。

圖4 矩形R中點(diǎn)的稀疏性
若矩形R中有多于6個(gè)S中的點(diǎn),則由鴿舍原理易知至少有一個(gè)δ×2δ的小矩形中有2個(gè)以上S中的點(diǎn)。設(shè)u,v是這樣2個(gè)點(diǎn),它們位于同一小矩形中,則

因此d(u,v)≤5δ/6<δ 。這與δ的意義相矛盾。也就是說(shuō)矩形R中最多只有6個(gè)S中的點(diǎn)。圖4(b)是矩形R中含有S中的6個(gè)點(diǎn)的極端情形。由于這種稀疏性質(zhì),對(duì)于P1中任一點(diǎn)p,P2中最多只有6個(gè)點(diǎn)與它構(gòu)成最接近點(diǎn)對(duì)的候選者。因此,在分治法的合并步驟中,我們最多只需要檢查6×n/2=3n對(duì)候選者,而不是n2/4對(duì)候選者。這是否就意味著我們可以在O(n)時(shí)間內(nèi)完成分治法的合并步驟呢?現(xiàn)在還不能作出這個(gè)結(jié)論,因?yàn)槲覀冎恢缹?duì)于P1中每個(gè)S1中的點(diǎn)p最多只需要檢查P2中的6個(gè)點(diǎn),但是我們并不確切地知道要檢查哪6個(gè)點(diǎn)。為了解決這個(gè)問(wèn)題,我們可以將p和P2中所有S2的點(diǎn)投影到垂直線(xiàn)l上。由于能與p點(diǎn)一起構(gòu)成最接近點(diǎn)對(duì)候選者的S2中點(diǎn)一定在矩形R中,所以它們?cè)谥本€(xiàn)l上的投影點(diǎn)距p在l上投影點(diǎn)的距離小于δ。由上面的分析可知,這種投影點(diǎn)最多只有6個(gè)。因此,若將P1和P2中所有S的點(diǎn)按其y坐標(biāo)排好序,則對(duì)P1中所有點(diǎn)p,對(duì)排好序的點(diǎn)列作一次掃描,就可以找出所有最接近點(diǎn)對(duì)的候選者,對(duì)P1中每一點(diǎn)最多只要檢查P2中排好序的相繼6個(gè)點(diǎn)。
至此,我們可以給出用分治法求二維最接近點(diǎn)對(duì)的算法CPAIR2如下:
function CPAIR2(S);
begin
if |S|=2 then δ:=S中這2點(diǎn)的距離
else if |S|=0
then δ:=∞
else begin
1. m:=S中各點(diǎn)x坐標(biāo)值的中位數(shù);
構(gòu)造S1和S2,使S1={p∈S|px≤m}和S2={p∈S|px>m}
2. δ1:=CPAIR2(S1);δ2:=CPAIR2(S2);
3. δm:=min(δ1,δ2);
4. 設(shè)P1是S1中距垂直分割線(xiàn)l的距離在δm之內(nèi)的所有點(diǎn)組成的集合,
P2是S2中距分割線(xiàn)l的距離在δm之內(nèi)所有點(diǎn)組成的集合。將P1和
P2中的點(diǎn)依其y坐標(biāo)值從小到大排序,并設(shè)P1*和P2*是相應(yīng)的已排
好序的點(diǎn)列;
5. 通過(guò)掃描P1*以及對(duì)于P1*中每個(gè)點(diǎn)檢查P2*中與其距離在δm之內(nèi)的
所有點(diǎn)(最多6個(gè))可以完成合并。當(dāng)P1*中的掃描指針逐次向上移動(dòng)
時(shí),P2*中的掃描指針可在寬為2δm的一個(gè)區(qū)間內(nèi)移動(dòng)。設(shè)δl是按
這種掃描方式找到的點(diǎn)對(duì)間的最小距離;
6. δ=min(δm,δl);
end;
return(δ);
end;
下面我們來(lái)分析一下算法CPAIR2的計(jì)算復(fù)雜性。設(shè)對(duì)于n個(gè)點(diǎn)的平面點(diǎn)集S,算法耗時(shí)T(n)。算法的第1步和第5步用了O(n)時(shí)間,第3步和第6步用了常數(shù)時(shí)間,第2步用了2T(n/2)時(shí)間。若在每次執(zhí)行第4步時(shí)進(jìn)行排序,則在最壞情況下第4步要用O(nlogn)時(shí)間。這不符合我們的要求。因此,在這里我們要作一個(gè)技術(shù)上的處理。我們采用設(shè)計(jì)算法時(shí)常用的預(yù)排序技術(shù),即在使用分治法之前,預(yù)先將S中n個(gè)點(diǎn)依其y坐標(biāo)值排好序,設(shè)排好序的點(diǎn)列為P*。在執(zhí)行分治法的第4步時(shí),只要對(duì)P*作一次線(xiàn)性?huà)呙瑁纯沙槿〕鑫覀兯枰呐藕眯虻狞c(diǎn)列P1*和P2*。然后,在第5步中再對(duì)P1*作一次線(xiàn)性?huà)呙瑁纯汕蟮?em>δl。因此,第4步和第5步的兩遍掃描合在一起只要用O(n)時(shí)間。這樣一來(lái),經(jīng)過(guò)預(yù)排序處理后的算法CPAIR2所需的計(jì)算時(shí)間T(n)滿(mǎn)足遞歸方程:

顯而易見(jiàn)T(n)=O(nlogn),預(yù)排序所需的計(jì)算時(shí)間為O(n1ogn)。因此,整個(gè)算法所需的計(jì)算時(shí)間為O(nlogn)。在漸近的意義下,此算法已是最優(yōu)的了。
晚上研究了半天,又參考了類(lèi)似的代碼,終于搞出來(lái)了:
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>

const int MAXN = 100001;
const double eps = 1e-6;

struct point
{
int index;
double x,y;
}a[MAXN],b[MAXN],c[MAXN];


inline double min(const double p,const double q)
{
return p>q ? q:p;
}

inline double distance(const point &p,const point &q)
{
return sqrt((p.x-q.x)*(p.x-q.x)+(p.y-q.y)*(p.y-q.y));
}

int cmpx(const void *p,const void *q)
{
double t = ((point *)p)->x - ((point *)q)->x;
if(t>eps) return 1;
else if(fabs(t)<=eps) return 0;
else return -1;
}

int cmpy(const void *p,const void *q)
{
double t = ((point *)p)->y - ((point *)q)->y;
if(t>eps) return 1;
else if(fabs(t)<=eps) return 0;
else return -1;
}

void merge(point p[],point q[],int s,int m,int t)
{
int i=s,j=m+1,k=s;

while(i<=m && j<=t)
{
if(q[i].y>q[j].y) p[k++]=q[j++];
else p[k++]=q[i++];
}
while(i<=m) p[k++]=q[i++];
while(j<=t) p[k++]=q[j++];
}

double closest(point a[],point b[],point c[],int p,int q)
{
if(q==p+1) return distance(a[p],a[q]);

if(q==p+2)
{
double d1=distance(a[p],a[q]);
double d2=distance(a[p],a[p+1]);
double d3=distance(a[p+1],a[q]);
if(d1<d2 && d1<d3) return d1;
else if(d2<d3) return d2;
else return d3;
}
int i,j,k,m=(p+q)>>1;
double d1,d2,dm;

for(i=j=p,k=m+1;i<=q;i++)
{
if(b[i].index<=m) c[j++]=b[i];
else c[k++]=b[i];
}
d1=closest(a,c,b,p,m),d2=closest(a,c,b,m+1,q);
dm=min(d1,d2);
merge(b,c,p,m,q);
for(i=k=p;i<=q;i++)
if(fabs(b[i].x-b[m].x)<dm)
c[k++]=b[i];
for(i=p;i<k;i++)

for(j=i+1;j<k && (c[j].y-c[i].y)<dm;j++)
{
double t=distance(c[i],c[j]);
if(t<dm) dm=t;
}
return dm;
}

int main()
{
int i,n;

while(scanf("%d",&n),n)
{
for(i=0;i<n;i++) scanf("%lf %lf",&a[i].x,&a[i].y);
qsort(a,n,sizeof(a[0]),cmpx);
for(i=0;i<n;i++) a[i].index=i;
memcpy(b,a,n*sizeof(a[0]));
qsort(b,n,sizeof(b[0]),cmpy);
printf("%.2lf\n",closest(a,b,c,0,n-1));
}
return 0;
}