算法导论-分治法-最近点对-HDOJ1007

HDOJ1007的原题目是求出在不同时套中两个玩具的前提下，圆圈的最大半径。问题翻译过来就是求解最近点对的问题，这个问题是经典的分治法问题。

参考博客：http://www.cnblogs.com/peng-come-on/archive/2012/01/18/2325163.html

毫无疑问，通过暴力手段列举所有的点对并计算这些点对的距离，找出最小的一组，可以得到最后的结果。但是，这道题的数据规模非常大，所以，这种传统的方法肯定行不通。我试过，hdoj是超时的。

我对原博客的第3步和第4步有不同的看法。

原博客步骤：

1.读取数据，并将点按横坐标升序排列。
2.以最中间的那个点为基准，将平面内的点分为左右两个部分。递归调用mindis(int,int)函数，分别求出左右两个部分的点集的最短距离，并取两者中的较小值，即为min.
3.显然，min并不一定是最短距离，因为还可能存在一种情况，即点对中的一个点位于左区域，另一个点位于右区域。所以，我们取点集中横坐标与分界线的距离小于min的点，存入p2[N]数组中。
4.对p2[N]数组中的点按纵坐标进行排序，计算p2[N]数组中的点对的距离，如果存在小于min 的情况，就取代min作为最近距离。
最后，min即为平面内点集的最近距离。

我的步骤是：

1.读取数据，并将点按横坐标升序排列。
2.以最中间的那个点为基准，将平面内的点分为左右两个部分。递归调用mindis(int,int)函数，分别求出左右两个部分的点集的最短距离，并取两者中的较小值，即为min.
3.显然，min并不一定是最短距离，因为还可能存在一种情况，即点对中的一个点位于左区域，另一个点位于右区域。所以，我们取点集中横坐标与分界线的距离小于min的点，根据在分界线的左右，分别存入pxSmall[N]数组和pxLarge[N]数组中。

4.因为唯一可能的情况是一个点在pxSmall[N]，另一个点在pxLarge[N]中。只要遍历这个两个数组便可以了。

下面是在原博客代码上的修改：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>

using namespace std;

#define N 1000010

struct point 
{
 double x;
 double y;
}p1[N],pxSmall[N],pxLarge[N];

double dis ( point a , point b )
{
	return sqrt( pow (a.x-b.x,2) + pow ( a.y-b.y,2 ) );
}

double min ( double a , double b )
{
	return a<b?a:b;
}

bool cpx ( point a , point b )
{
	return a.x < b.x ;
}
bool cpy ( point a , point b )
{
	return a.y < b.y ;
}

double mindis (int l, int r)
{
	if( l + 1 == r )
		return dis ( p1[l] ,p1[r] );
	if( l + 2 == r )
		return min ( dis ( p1[l] , p1[l+1] ) , min ( dis ( p1[l+1] , p1[r] ) , dis ( p1[l] , p1[r] ) ) );
	else
	{
		int mid ,count1=0, count2=0;
		double mini;
		mid = ( l + r) >> 1 ;
		mini = min ( mindis ( l , mid ) , mindis ( mid+1 , r ) );
		for( int i = l ; i <= r ; i++ )
		{
			if ( fabs ( p1[i].x - p1[mid].x ) <= mini )
			{
				if (p1[i].x-p1[mid].x < 0)
					pxSmall[count1++]=p1[i];
				else
					pxLarge[count2++]=p1[i];
			}
				
		}
		//直接遍历两个数组
		for(int i=0;i<count1;i++)
		{
			for(int j=0;j<count2;j++)
			{
				double temp = dis(pxSmall[i], pxLarge[j]); 
				if(temp<mini)
					mini=temp;
			}
		}
		/*
		sort ( p2 , p2+count , cpy );
		for ( int i=0 ; i < count ; i++ )
		{
			for ( int j = i+1; j < count ;j++)
			{
				if ( p2[j].y-p2[i].y>=mini)
					break;
				else if(dis (p2[j],p2[i])<mini)
					mini=dis(p2[j],p2[i]);
			}
		}
		*/
		return mini;
	}
}
  
  

int main()
{
 //freopen("input.txt","r",stdin);
	int n ;
	double dia ;
	while(scanf("%d",&n)==1&&n)
	{
		for(int i=0;i<n;i++)
		scanf("%lf%lf",&p1[i].x,&p1[i].y);
		sort ( p1 , p1 + n-1 , cpx );
		dia = mindis ( 0 , n-1 );
		printf("%.2f\n", dia / 2 );
	}
	return 0;
}

我在HDOJ上分别跑这两组代码，发现原博客runtime为1875ms，我的代码runtime为1093ms，比原博客的运行时间少了很多。因为少了排序这个步骤。当然，可能不同的数据集得到的运行时间不同，如果换一个数据集或数据集规模很大的时候原博客的方法可能更有优势。但是有一点，我的方法更好地符合分治法的要求，正确性是可以保证的。最近一直在学习《算法导论》，这本书真是神书，是算法中的算法，不仅仅告诉你算法，还会告诉你算法的证明。以后如果遇到分治法的问题，我会继续补充到这篇博客。

补充：参考博客：http://blog.youkuaiyun.com/junerfsoft/article/details/2975495

博客中提到了飞机调度的问题，说的很好。可惜因为每个点的坐标类型都是double类型的，不能直接找到那6个点，可以使用二分查找，但是也不是特别方便。我随便写了代码，时间上反而更长了。