1.10 HDU 1007 Quoit Design
http://acm.hdu.edu.cn/showproblem.php?pid=1007
首先说明题意：一个名为quiot的游戏，类似与在商场玩的用圆环套玩具的游戏，但设计要求保证一次不应套中两个玩具。求计算最大的圆环半径。
于是可以建模为平面最近点对的经典分治问题。
1 #include <stdio.h> 2 #include <math.h> 3 #include <algorithm> 4 using namespace std; 5 const int MAXN = 100100; 6 struct Point{ 7 double x, y; 8 }x[MAXN]; 9 int y[MAXN], z[MAXN]; 10 bool cmpx(const Point& a, const Point& b) 11 { return a.x < b.x; } 12 bool cmpy(const int& a, const int &b) 13 { return x[a].y < x[b].y; } 14 inline double dis(Point a, Point b) 15 { return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y); } 16 17 void Merge(int z[], int y[], int l, int r, int m) 18 { 19 int i = l, j = m+1, k = l; 20 while(i <= m && j <= r) 21 if(x[z[i]].y <= x[z[j]].y) y[k++] = z[i++]; 22 else y[k++] = z[j++]; 23 while(i <= m) y[k++] = z[i++]; 24 while(j <= r) y[k++] = z[j++]; 25 } 26 27 double close(int y[], int z[], int l, int r) 28 { 29 if(l+1 == r) 30 return dis(x[l], x[r]); 31 if(l+2 == r) 32 return min(dis(x[l], x[l+1]), min(dis(x[l+1], x[r]), dis(x[l], x[r]))); 33 int m = l+r>>1; 34 for(int i=l, j=m+1, k=l ; k<=r ; k++) 35 if(y[k] > m) z[j++] = y[k]; 36 else z[i++] = y[k]; 37 double res = min(close(z, y, l, m), close(z, y, m+1, r)); 38 Merge(z, y, l, r, m); 39 int rt = l; 40 for(int i=l ; i<=r ; i++) 41 if(fabs(x[y[i]].x - x[y[m]].x) < res) 42 z[rt++] = y[i]; 43 for(int i=l ; i<rt ; i++) 44 for(int j=i+1 ; j<rt ; j++) 45 { 46 if(x[z[j]].y - x[z[i]].y >= res) break; 47 res = min(res, dis(x[z[i]], x[z[j]])); 48 } 49 return res; 50 } 51 52 int main() 53 { 54 // freopen("input", "r", stdin); 55 int n; 56 while(scanf("%d", &n), n) 57 { 58 for(int i=0 ; i<n ; i++) 59 scanf("%lf %lf", &x[i].x, &x[i].y), y[i] = i; 60 sort(x, x+n, cmpx); 61 sort(y, y+n, cmpy); 62 printf("%.2lf\n", sqrt(close(y, z, 0, n-1))/2); 63 } 64 }
1. X，Y，Z数组的功能：X数组存放了所有点记录，按照x坐标排序，在每一个分治函数中都不改变；Y数组用于开始也存放了所有点记录，按照y坐标排序，Z是临时数组。这里用了一个很巧妙的方法使空间复杂度降为O（n），即让YZ轮换。每次分治时，我们要保证传入的Y数组在相应[l, r]范围内是不变的。我们可以把[l, m]部分和[m+1, r]部分的点分别按照y升序放入Z数组相应的[l, r]，这样Y数组[l, r]部分就没有用了，可以让Y数组起到下一次分治临时数组的功能（代码反映在：close(z, y, l, m), close(z, y, m+1, r)）。但是当分治完需要回复Y数组（记得一开始我们说需要保证Y数组的不变性），而根据递归定义，传入的Z的[l, m]和[m+1, r]部分是不变的，故可以归并恢复Y数组。
2. 所谓d*2d内最多有六个点体现在那里？我们要用分治构造出整体O（nlogn）的算法，就要求符合T(n) = 2T(n/2) + O(n)的递归式，即close函数内只能进行O(n)的算法。34-36行把Y数组按照m分割为两个部分放入Z数组，38行Merge归并构造Y，39-42行将距离m两边d的所有点存入Z数组，均是一个for循环从l到r，故符合O(n)复杂度。43-48行有两个for循环，其功能是在候选点集中一一计算距离，正是这里体现了鸽笼定理证出的“六个点”。内部的break可以保证第二个循环至多执行3次。所以也是O（n）的复杂度。
3. 为什么不需要l == r的判断？由于我们单独处理的l+2 == r和l+1 == r，这样对r-l >= 3执行m = l+r>>1不可能出现m==l或m==r。
4. 还有一个需要注意的地方：我们在close函数中对Y和Z数组都限定了范围：l-r，即使Z是临时数组，我们可以利用的也只有[l, r]的范围。这样才有了算法的正确性。
5. 令我费解的一点是，如果把第60行的sort排序换成qsort，就会WA（改61行的不会WA）。。。这让我这个C爱好者感到好桑心，每次码题都强迫征似的只包含C的头文件，结果WA了无数次。。。这里贴上qsort的代码，可能是比较函数错了吧。。。求大神指出为何会错。
1 #include <stdio.h> 2 #include <string.h> 3 #include <stdlib.h> 4 #include <math.h> 5 const int MAXN = 100100; 6 struct Point{ 7 double x, y; 8 }x[MAXN]; 9 int y[MAXN], z[MAXN]; 10 int cmpx(const void* a, const void* b) 11 { return ((Point*)a)->x - ((Point*)b)->x; } 12 int cmpy(const void* a, const void* b) 13 { return x[(*(int*)a)].y - x[(*(int*)b)].y; } 14 inline double min(double a, double b) 15 { return a < b ? a : b; } 16 inline double dis(Point a, Point b) 17 { return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y); } 18 void Merge(int z[], int y[], int l, int r, int m) 19 { 20 int i = l, j = m+1, k = l; 21 while(i <= m && j <= r) 22 if(x[z[i]].y <= x[z[j]].y) y[k++] = z[i++]; 23 else y[k++] = z[j++]; 24 while(i <= m) y[k++] = z[i++]; 25 while(j <= r) y[k++] = z[j++]; 26 } 27 28 double close(int y[], int z[], int l, int r) 29 { 30 if(l+1 == r) 31 return dis(x[l], x[r]); 32 if(l+2 == r) 33 return min(dis(x[l], x[l+1]), min(dis(x[l+1], x[r]), dis(x[l], x[r]))); 34 int m = l+r>>1; 35 for(int i=l, j=m+1, k=l ; k<=r ; k++) 36 if(y[k] > m) z[j++] = y[k]; 37 else z[i++] = y[k]; 38 double res = min(close(z, y, l, m), close(z, y, m+1, r)); 39 Merge(z, y, l, r, m); 40 int rt = l; 41 for(int i=l ; i<=r ; i++) 42 if(fabs(x[y[i]].x - x[y[m]].x) < res) 43 z[rt++] = y[i]; 44 for(int i=l ; i<rt ; i++) 45 for(int j=i+1 ; j<rt ; j++) 46 { 47 if(x[z[j]].y - x[z[i]].y >= res) break; 48 res = min(res, dis(x[z[i]], x[z[j]])); 49 } 50 return res; 51 } 52 53 int main() 54 { 55 // freopen("input", "r", stdin); 56 int n; 57 while(scanf("%d", &n), n) 58 { 59 for(int i=0 ; i<n ; i++) 60 scanf("%lf %lf", &x[i].x, &x[i].y), y[i] = i; 61 qsort(x, n, sizeof(Point), cmpx); 62 qsort(y, n, sizeof(int), cmpy); 63 printf("%.2lf\n", sqrt(close(y, z, 0, n-1))/2); 64 } 65 }
当然我还看了网上的做法：我看到网上的好多方法并没有用Merge操作，代码比较短，于是自己照着写了一遍
1 #include <stdio.h> 2 #include <string.h> 3 #include <stdlib.h> 4 #include <math.h> 5 #include <algorithm> 6 using namespace std; 7 #define dis(a, b) (((a->x-b->x)*(a->x-b->x) + (a->y-b->y)*(a->y-b->y))) 8 const int MAXN = 100100; 9 struct Point{ 10 double x, y; 11 }p[MAXN], *px[MAXN], *py[MAXN]; 12 bool cmpx( Point *p1, Point *p2 ) 13 { return p1->x < p2->x; } 14 bool cmpy( Point *p1, Point *p2 ) 15 { return p1->y < p2->y; } 16 17 double close(int l, int r) 18 { 19 if(l+1 == r) 20 return dis(px[l], px[r]); 21 if(l+2 == r) 22 return min(dis(px[l], px[l+1]), min(dis(px[l+1], px[r]), dis(px[l], px[r]))); 23 int m = l+r>>1, rt = 0; 24 double d = min(close(l, m), close(m+1, r)); 25 for(int i=l ; i<=r ; i++) 26 if(px[i]->x >= px[m]->x - d && px[i]->x <= px[m]->x + d) 27 py[rt++] = px[i]; 28 sort(py, py+rt, cmpy); 29 for(int i=0 ; i<rt ; i++) 30 for(int j=i+1 ; j<rt ; j++) 31 { 32 if(py[j]->y - py[i]->y >= d) break; 33 d = min(d, dis(py[i], py[j])); 34 } 35 return d; 36 } 37 38 int main() 39 { 40 // freopen("input", "r", stdin); 41 int n; 42 while(scanf("%d", &n), n) 43 { 44 for(int i=0 ; i<n ; i++) 45 scanf("%lf %lf", &p[i].x, &p[i].y), px[i] = &p[i]; 46 sort(px, px+n, cmpx); 47 printf("%.2lf\n", sqrt(close(0, n-1))/2.0); 48 } 49 }
从评测结果上看，比我的快几十毫秒吧，但从复杂度上来说，这是个更渣的算法。在第28行出close函数里用了一个排序函数，虽然表面上看只排了rt个元素，但最坏情况下rt会等于r-l，这样递归方程变成了T(n) = 2T(n/2) + O(nlogn)，根据主定理，总复杂度O(n) = nlognlogn。从复杂度上分析这种方式更糟糕，但hdu上结果跟nlogn的结果差不多，可能还是数据太小的缘故吧。
不知为何，把sort换为qsort一样会WA，而且是在我仔细看了qsort的用法后。。。看来开学需要向人好好请教下qsort的用法了。
