1.10 HDU 1007 Quoit Design

http://acm.hdu.edu.cn/showproblem.php?pid=1007

首先说明题意:一个名为quiot的游戏,类似与在商场玩的用圆环套玩具的游戏,但设计要求保证一次不应套中两个玩具。求计算最大的圆环半径。

于是可以建模为平面最近点对的经典分治问题。

1.10 HDU 1007 Quoit Design
 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 }
View Code

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.10 HDU 1007 Quoit Design
 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 }
View Code

 

当然我还看了网上的做法:我看到网上的好多方法并没有用Merge操作,代码比较短,于是自己照着写了一遍

1.10 HDU 1007 Quoit Design
 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 }
View Code

从评测结果上看,比我的快几十毫秒吧,但从复杂度上来说,这是个更渣的算法。在第28行出close函数里用了一个排序函数,虽然表面上看只排了rt个元素,但最坏情况下rt会等于r-l,这样递归方程变成了T(n) = 2T(n/2) + O(nlogn),根据主定理,总复杂度O(n) = nlognlogn。从复杂度上分析这种方式更糟糕,但hdu上结果跟nlogn的结果差不多,可能还是数据太小的缘故吧。

不知为何,把sort换为qsort一样会WA,而且是在我仔细看了qsort的用法后。。。看来开学需要向人好好请教下qsort的用法了。

1.10 HDU 1007 Quoit Design

上一篇:CSS3 Box-sizing


下一篇:[译] 第二十七天:Restify - 用Node.js创建正确的REST Web服务