问题描述:以(x,y)表示一个点,随机生成n个点,求出最近的一对点。
解法一:蛮力法
思路:用一个变量minD记录最近的距离,初始值为无穷大(int的最大表示范围),然后在n个点选取2个点(组合问题,一共有C n取2种可能),计算它们的距离,如果比minD小,更新minD的值,一直循环,复杂度为O(n*n)。
上代码(下述所有代码样例的头文件都一样,用codeblock编译器执行):
#include <iostream>
#include<random>
#include<ctime>
#include<string>
#include <limits>
#include<vector>
#include<math.h>
#include<algorithm>
using namespace std;
class point
{
public:
int x,y;
point(int a,int b){x=a,y=b;}
point(){x=y=0;}
void Copy(point aPoint)
{
x=aPoint.x;
y=aPoint.y;
}
};
void generateNPoints(int n,point * points)
{
//c++11均匀分布
uniform_int_distribution<unsigned> u(0,10000);//点的生成分为是0~10000,多了怕计算距离的时候溢出
//随机数引擎
default_random_engine e( clock() );
for(int i=0;i<n;++i)
{
points[i].x = u(e);
points[i].y = u(e);
}
}
//打印顶点数组
void print(int n,point * points)
{
for(int i=0;i<n;++i)
{
cout<<'('<<points[i].x<<','<<points[i].y<<')'<<' ';
if((i+1)%5==0)
cout<<endl;
}
cout<<endl;
}
//平方
int pow2(int x)
{
return x*x;
}
//蛮力法求解最近点对问题
int force(int n,point * points)
{
point point1,point2;
int i,j,d,minD;
d=minD = (numeric_limits<int>::max)();
//遍历每一个点
for(i=0;i<n;++i)
{
//每一个点都跟它后面的点比较
for(j=i+1;j<n;++j)
{
//计算距离的平方
d = pow2(points[i].x-points[j].x) + pow2(points[i].y-points[j].y);
//如果距离比minD小,更新minD
if(d<minD)
{
minD = d;
//point1,point2用来存放点对
point1.Copy(points[i]);
point2.Copy(points[j]);
}
}
}
cout<<"蛮力法求得最近点对为("<<point1.x<<','<<point1.y<<"),("<<point2.x<<','<<point2.y<<")。"<<endl;
return sqrt(minD);
}
int main()
{
int n=100,t1=1,t2,i;
clock_t startTime;
clock_t endTime ;
int time[20];
double aveTime;
while(t1--)
{
n*=10;
t2= 20;
i = 0;
aveTime = 0.0;
while(t2--)
{
minD =(numeric_limits<int>::max)();
point * points1;
points1 = new point[n];
generateNPoints(n,points1);
startTime = clock();
cout<<"蛮力法求得最短距离为:"<<force(n,points1)<<endl;
endTime = clock();
time[i++] = double(endTime-startTime)*1000/CLOCKS_PER_SEC;
delete [] points1;
}
cout<<"scale = "<<n<<endl;
for(i=0,aveTime =0;i<40;++i)
{
aveTime += time[i];
}
aveTime *= (1/20.0);
cout<<"使用蛮力法平均运行时间为"<<aveTime<<"ms"<<endl;
}
return 0;
}
运行结果:
本人还测试过,当规模为100,000时,运行时间为38秒左右。
解法二:分治法
相信大家已经看过为什么分治法可以解决这个问题,我就不解释了,因为我也是看到网上的==,下面只讲重点。
显然重点是分治以后的合并,请看下图:
一种思路是先找出处于两条绿色虚线中间的点,然后对y轴排序,根据抽屉原理,每个点最多只需对后面的7个点遍历(具体理由见算法导论610页)。性能分析:对点按y轴排序最少需要O(n*lgn)步骤,遍历需要7n,故算法复杂度是O(n* (lgn)^2).
核心代码:
//得到两点距离的平方
int getD2(point p1,point p2)
{
return pow2(p1.x-p2.x) + pow2(p1.y-p2.y);
}
//按x值比较
bool cmpX(const point &p1,const point &p2)
{
return p1.x < p2.x;
}
//按y值比较
bool cmpY(const point &p1,const point &p2)
{
return p1.y < p2.y;
}
int minD ;
point point1, point2;
//递归调用分治函数
int Devide(point * points1,int be,int en)
{
//当只有一个点时返回无穷大
if(be >= en)
{
return (numeric_limits<int>::max)();
}
//当有两个点时直接返回他们的距离
else if(en - be == 1)
{
int d=sqrt( getD2(points1[be],points1[en]) );
if(d<minD)
{
minD = d;
point1.Copy(points1[be]);
point2.Copy(points1[en]);
}
return d;
}
//否则当前点集按下标值划分成两部分
else//be - en>1
{
int d1,d2,d3,mid,d;
vector<point> s;
mid = (en + be)/2;
//求得左边的最短距离
d1 = Devide(points1,be,mid);
//求得右边的最短距离
d2 = Devide(points1,mid+1,en);
d=d3= (numeric_limits<int>::max)();
d = d1>d2?d2:d1;
minD = d<minD?d:minD;
//从分割线向两边寻找虚线内的点集,符合就放到向量s里去
for(int i = mid;i >= be && points1[mid].x-points1[i].x <= d ; --i)
{
s.push_back(points1[i]);
}
for(int i = mid+1;i <= en && points1[i].x - points1[mid].x <= d ; ++i)
{
s.push_back(points1[i]);
}
//对符合要求的点集按y值排序
sort(s.begin(),s.end(),cmpY);
int l = s.size();
//对于s里的每一个点
for(int i = 0; i < l ; ++i)
{
//最多遍历7次
int limit = i+8 > l ? l :i+8;
for(int j = i+1; j<limit && s[j].y - s[i].y < d ; ++j)
{
//计算两点距离
d3 = sqrt( getD2(s[j],s[i]) );
//如果比minD小,更新点对和minD
if(d3 < minD)
{
minD = d3;
point1.Copy(s[i]);
point2.Copy(s[j]);
}
}
}
return minD;
}
}
//分治法求解最近点对问题,points1是按x值排序后的顶点数组
int devide(int n,point * points1)
{
int d;
d =Devide(points1,0,n-1);
cout<<"分治法求得最近点对为("<<point1.x<<','<<point1.y<<"),("<<point2.x<<','<<point2.y<<")。"<<endl;
return d;
}
运行结果:
分治法的优化:
通过上面的分析我们知道,耗时的操作在于每次递归的时候都要对虚线内的点集按y值排序,那么我们是否可以在递归之前就排序?我们同时传一个按x值和y值排序的数组过来,每次递归的时候遍历一次按y值排序的数组,寻找虚线里面的点集,其余跟上面一致。注意每次递归都要更新按y值排序的数组,该数组只需要包含当前范围内的点。
上代码:
int minD ;
point point1,point2;
int Devide(int n,point * points1,point * points2,int be,int en)
{
if(be >= en)
{
return (numeric_limits<int>::max)();
}
else if(en - be == 1)
{
int d=sqrt( getD2(points1[be],points1[en]) );
if(d<minD)
{
minD = d;
point1.Copy(points1[be]);
point2.Copy(points1[en]);
}
return d;
}
else//be - en>1
{
int d1,d2,d3,mid,d,leftNum=0,i,j,k;
vector<point> s;
mid = (en + be)/2;
//更新按y值排序的数组
for(i =0;i<n;++i)
{
if(points1[mid].x-points2[i].x >=0)
++leftNum;
}
point * pl,*pr;
pl = new point [leftNum];
pr = new point [n - leftNum];
for(i=j=k=0;i<n;++i)
{
if(points1[mid].x-points2[i].x >=0)
pl[j++].Copy(points2[i]);
else
pr[k++].Copy(points2[i]);
}
//递归调用,一分为二
d1 = Devide(j,points1,pl,be,mid);
d2 = Devide(k,points1,pr,mid+1,en);
delete [] pl;
delete [] pr;
d=d3= (numeric_limits<int>::max)();
d = d1>d2?d2:d1;
minD = d<minD?d:minD;
//从分割线向两边寻找虚线内的点集,符合就放到向量s里去
for(int i = 0;i <n ; ++i)
{
if(abs( points2[i].x - points1[mid].x) <= d)
s.push_back(points2[i]);
}
int l = s.size();
for(int i = 0; i < l ; ++i)
{
int limit = i+8 > l ? l :i+8;
for(int j = i+1; j<limit && s[j].y - s[i].y < d ; ++j)
{
d3 = sqrt( getD2(s[j],s[i]) );
if(d3 < minD)
{
minD = d3;
point1.Copy(s[i]);
point2.Copy(s[j]);
}
}
}
return minD;
}
}
//分治法求解最近点对问题,points1按x值排序,points2按y值排序
int devide(int n,point * points1,point * points2)
{
int d;
d =Devide(n,points1,points2,0,n-1);
cout<<"分治法求得最近点对为("<<point1.x<<','<<point1.y<<"),("<<point2.x<<','<<point2.y<<")。"<<endl;
return d;
}
运行结果:
此外还有一个可以优化的地方,就是递归函数的返回条件可以再增加一个:只有三个点的时候。
//写在递归函数Devide的返回条件里
else if(en - be == 2)
{
int d1,d2,d3,mid,d;
mid = (en + be)/2;
d1 = getD2(points1[be],points1[mid]);
d2 = getD2(points1[be],points1[en]);
d3 = getD2(points1[mid],points1[en]);
if(d1 <= d3 && d1 <= d2)
{
d =sqrt(d1);
if(d < minD)
{
minD = d;
point1.Copy(points1[be]);
point2.Copy(points1[mid]);
}
return d;
}
else if(d2 <= d3 && d2 <= d1)
{
d =sqrt(d2);
if(d < minD)
{
minD = d;
point1.Copy(points1[be]);
point2.Copy(points1[en]);
}
return d;
}
else
{
d =sqrt(d3);
if(d < minD)
{
minD = d;
point1.Copy(points1[en]);
point2.Copy(points1[mid]);
}
return d;
}
}
运行结果:
我觉得200ms左右应该是比较理想的状态了。
在规模为1000,000下,运行时间为一秒八左右。