快排、堆排对于随机数据的效率测试

放假期间无所事事，决定回顾一下快速排序和堆排序。
自从转C后，快排都用qsort解决，去年参加CCF培训的时候有些选手说stdlib里面的qsort效率低，自己写更好。
堆排序实际并不常用，除非用于一些变动频繁的数据，例如dijstra算法中的最小堆优化。
堆排和快排的时间复杂度都是O(nlogn)，但这并不代表它们的实际效率是相同的，于是就有了下面的测试。

先介绍下环境：CPU~Dual core Intel Core i3 CPU M 370 (-HT-MCP-) clocked at 933.000 Mhz Kernel~3.0.0-17-generic x86_64 Up~4:08 Mem~674.2/3760.7MB HDD~320.1GB(44.1% used) Procs~181 Client~Shell inxi~1.7.7
首先，随机生成10000000个数字，存到rand.txt

然后，用下面的代码测试：

C语言: 高亮代码由发芽网提供

001 #include
002 #include
003 #include
004 #include
005 #define MAXN 10000000
006
007 long num[MAXN];
008 long bak[MAXN];
009 long heapsize;
010
011 long lc(long k)
012 {
013     return 2 * k + 1;
014 }
015
016 long rc(long k)
017 {
018     return 2 * k + 2;
019 }
020
021 long p(long k)
022 {
023     return (k - 1) >> 1;
024 }
025
026 void read()
027 {
028     long i;
029     for(i = 0;i < MAXN;i ++)
030         scanf("%ld",&num[i]);
031 }
032
033 void pt()
034 {
035     long i;
036     for(i = 0;i < MAXN;i ++)
037         printf("%ld\n",num[i]);
038 }
039
040 int cmp(const void *a,const void *b)
041 {
042     return *(long *)a > *(long *)b;
043 }
044
045 void quicksort(long l,long r)
046 {
047     if(l < r)
048     {
049     long mid,key,tmp,i,j;
050     mid = (l + r) >> 1;
051     key = num[mid];num[mid] = num[r];num[r] = key;
052     for(i = j = l;i < r;i ++)
053     {
054         if(num[i] < key)
055         {
056             tmp = num[i];
057             num[i] = num[j];
058             num[j] = tmp;
059             j ++;
060         }
061     }
062     num[r] = num[j];
063     num[j] = key;
064     quicksort(l,j - 1);
065     quicksort(j + 1,r);
066     }
067 }
068
069 void quicksort_rand(long l,long r)
070 {
071     if(l < r)
072     {
073     long mid,key,tmp,i,j;
074     mid = rand() % (r - l) + l;
075     key = num[mid];num[mid] = num[r];num[r] = key;
076     for(i = j = l;i < r;i ++)
077     {
078         if(num[i] < key)
079         {
080             tmp = num[i];
081             num[i] = num[j];
082             num[j] = tmp;
083             j ++;
084         }
085     }
086     num[r] = num[j];
087     num[j] = key;
088     quicksort(l,j - 1);
089     quicksort(j + 1,r);
090     }
091 }
092
093 void heapfy(int k)
094 {
095     long min,tmp;
096     min = k;
097     if(lc(k) < heapsize && num[lc(k)] < num[min])    min = lc(k);
098     if(rc(k) < heapsize && num[rc(k)] < num[min])    min = rc(k);
099     if(min != k)
100     {
101         tmp = num[min];
102         num[min] = num[k];
103         num[k] = tmp;
104         heapfy(min);
105     }
106 }
107
108 long extract_min()
109 {
110     long min = num[0];
111     num[0] = num[--heapsize];
112     heapfy(0);
113     return min;
114 }
115
116 void heapsort()
117 {
118     long i;
119     for(i = MAXN / 2;i >= 0;i --)
120         heapfy(i);
121     while(heapsize >= 1)
122     {
123         i = extract_min();
124         //printf("%ld\n",i);
125     }
126 }
127
128 void create()
129 {
130     FILE *out;
131     out = fopen("rand.txt","w");
132     long i;
133     srand(time(NULL));
134     for(i = 0;i < MAXN;i ++)
135         fprintf(out,"%d\n",rand());
136     exit(0);
137 }
138
139 int main()
140 {
141     clock_t start,end;
142     //create();
143     read();
144     memcpy(bak,num,sizeof(num));
145     start = clock();
146     qsort(num,MAXN,sizeof(long),cmp);
147     end = clock();
148     //pt();
149     fprintf(stderr,"Biult-in qsort:%f\n",(double)(end - start) / CLOCKS_PER_SEC);
150
151     memcpy(num,bak,sizeof(num));
152     start = clock();
153     quicksort(0,MAXN - 1);
154     end = clock();
155     fprintf(stderr,"My QuickSort:%f\n",(double)(end - start) / CLOCKS_PER_SEC);
156     //pt();
157    
158     memcpy(num,bak,sizeof(num));
159     start = clock();
160     quicksort_rand(0,MAXN - 1);
161     end = clock();
162     fprintf(stderr,"My QuickSort(Rand):%f\n",(double)(end - start) / CLOCKS_PER_SEC);
163     //pt();
164    
165     memcpy(num,bak,sizeof(num));
166     start = clock();
167     heapsize = MAXN;
168     heapsort();
169     end = clock();
170     fprintf(stderr,"HeapSort:%f\n",(double)(end - start) / CLOCKS_PER_SEC);
171     return 0;
172 }

当MAXN定义为1000000时，输出如下：

Biult-in qsort:0.250000
My QuickSort:0.240000
My QuickSort(Rand):0.230000
HeapSort:0.830000

当MAXN达到10000000时：
Biult-in qsort:3.010000
My QuickSort:2.640000
My QuickSort(Rand):2.640000
HeapSort:13.430000

整体上，可以看出，对于随机数据，单次排序，快排的表现明显优于推排序。个人推测原因是堆排序存在建堆、弹出小元素和维护堆三个过程，因此耗时会多一些。

对于快速排序，qsort使用方便，效率略低。自己写的随机化快排速度最快。不过三者差别并不明显，本人还是推荐用stdlib内置的qsort——简单，安全。

输入输出效率分析

请先看这篇博文，注意字符串的读入方式：http://jim-think.blogspot.com/2011/10/contact.html
在USACO平台上，上面博文中的代码运行最长运行时间～0.4s

在同学的启发下，本人修改了读入部分，换成了如下代码：

C语言: 高亮代码由发芽网提供

1 char c;
2 while(!feof(stdin))
3 {
4     c = getchar();
5     if(c == '0')    str[len ++] = 0;
6     else if(c == '1')    str[len ++] = 1;
7 }

重新提交，结果如下：

   Test 1: TEST OK [0.000 secs, 11992 KB]
   Test 2: TEST OK [0.000 secs, 11992 KB]
   Test 3: TEST OK [0.000 secs, 11992 KB]
   Test 4: TEST OK [0.000 secs, 11992 KB]
   Test 5: TEST OK [0.022 secs, 11992 KB]
   Test 6: TEST OK [0.022 secs, 11992 KB]
   Test 7: TEST OK [0.022 secs, 11992 KB]

程序运行时间大大缩短！
可见，strcat的效率十分低下。
In conclusion，选择合适的输入输出方式避免时间浪费。

USACO Snail Trail

一开始想复杂了，甚至拿来建图。其实直接用DFS模拟蜗牛的路径就可以了。貌似DFS效率并不是很低，最长耗时不过0.02s。

C语言: 高亮代码由发芽网提供

01 #include <stdio.h>
02 #include <stdlib.h>
03 #include <string.h>
04
05 int N;
06 int map[125][125];
07 int dx[4] = {0,0,-1,1};
08 int dy[4] = {1,-1,0,0};
09 int max;
10
11 void init()
12 {
13     int B,i;
14     char c;
15     scanf("%d%d\n",&N,&B);
16     while(B --)
17     {
18         scanf("%c%d\n",&c,&i);
19         map[i][(int)(c - 'A' + 1)] = 1;
20     }
21     for(i = 0;i <= N + 1;i ++)  //borders as blocks
22         map[0][i] = map[N + 1][i] = map[i][0] = map[i][N + 1] = 1;
23 }
24
25 int visit[125][125] = {};
26
27 void dfs(int i,int j,int step,int d)
28 {
29     if(step > max)    max = step;
30     if(map[i + dx[d]][j + dy[d]] == 1)  //block ahead, change direction
31     {
32         if(d <= 1)
33         {
34             int k;
35             for(k = 2;k <= 3;k ++)
36             if(map[i + dx[k]][j + dy[k]] == 0 && visit[i + dx[k]][j + dy[k]] == 0)
37             {
38                 visit[i + dx[k]][j + dy[k]] = 1;
39                 dfs(i + dx[k],j + dy[k],step + 1,k);
40                 visit[i + dx[k]][j + dy[k]] = 0;
41             }
42         }
43         if(d >= 2)
44         {
45             int k;
46             for(k = 0;k <= 1;k ++)
47             if(map[i + dx[k]][j + dy[k]] == 0 && visit[i + dx[k]][j + dy[k]] == 0)
48             {
49                 visit[i + dx[k]][j + dy[k]] = 1;
50                 dfs(i + dx[k],j + dy[k],step + 1,k);
51                 visit[i + dx[k]][j + dy[k]] = 0;
52             }
53         }
54     }
55     else
56     {
57         if(visit[i + dx[d]][j + dy[d]] == 0)   //go ahead
58         {
59             visit[i + dx[d]][j + dy[d]] = 1;
60             dfs(i + dx[d],j + dy[d],step + 1,d);
61             visit[i + dx[d]][j + dy[d]] = 0;
62         }
63     }
64 }
65
66 int main()
67 {
68     freopen("snail.in","r",stdin);
69     freopen("snail.out","w",stdout);
70     init();
71     memset(visit,0,sizeof(visit));
72     dfs(1,1,1,0);
73     memset(visit,0,sizeof(visit));
74     dfs(1,1,1,3);
75     printf("%d\n",max);
76     return 0;
77 }

USACO Music Theme

最近几天生病在家，因此做题速度略有上升。

这道题有两种思路：
1、枚举最长theme长度，然后寻找是否存在。可以用二分提升速度
2、动态规划
先说第一种吧。先处理出一个diff数组，diff[i] = note[i + 1] - note[i]，那么如果存在长度为n的，中间至少间隔1个元素的diff子串，那么就有长度为n+1的theme。所以就有了下面的代码。

C语言: 高亮代码由发芽网提供

01 /*
02 LANG:C
03 ID:niat8181
04 PROB:theme
05 */
06
07 #include <stdio.h>
08 #include <stdlib.h>
09 #include <string.h>
10
11 int note[5050];
12 int diff[5050];
13 long diff_sum[5050];   //for qiuckness in comparison
14 int n;
15
16 int find(int len)
17 {
18     int i,j;
19     int t;
20     for(i = 0;i <= n - 2 * len;i ++)
21     {
22         t = diff_sum[i + len] - diff_sum[i];
23         for(j = i + len;j < n - len;j ++)
24             if(diff_sum[j + len] - diff_sum[j] == t && memcmp(diff + i,diff + j,len * sizeof(int)) == 0)        return 1;
25     }
26     return 0;
27 }
28
29 int main()
30 {
31     freopen("theme.in","r",stdin);
32     freopen("theme.out","w",stdout);
33     scanf("%d\n",&n);
34     int i;
35     for(i = 0;i < n;i ++)
36         scanf("%d",&note[i]);
37     for(i = 0;i < n - 1;i ++)
38         diff[i] = note[i + 1] - note[i];
39     diff_sum[0] = 0;
40     for(i = 0;i < n - 1;i ++)
41         diff_sum[i + 1] = diff_sum[i] + diff[i];
42    
43     int l,r,mid;
44     r = n / 2 - 1;
45     l = 3;
46     while(l < r)  //二分
47     {
48         mid = (l + r + 1) / 2;
49         if(find(mid))
50             l = mid;
51         else
52             r = mid - 1;
53     }
54     if(l >= 4)    printf("%d\n",l + 1);   //there is a theme
55     else printf("0\n");   //no theme exists
56     return 0;
57 }

现在说动态规划的。话说俺好久没写动态规划题目了，思路没有以前那么敏捷。这个题目可以说是最长公共子串（我觉得有必要翻一下算法导论额）的改版把。
用len[i][j]表示分别以note[i],note[j]开头的乐曲的最长相似长度。那么:
len[i][j] = len[i + 1][j + 1] + 1
当note[i] - note[j] == note[i + 1] - note[j + 1]（即可以连接成一个更长的theme）；
否则len[i][j] = 1(自身长度)

读者是否注意到，这个动态规划形式是二维的，5000*5000的int数组需要约50M空间。那么，就此放弃动态规划吗？

在状态转移方程中可以发现，len[i][j]仅和len[i+1][j+1]有关系。这样，一维数组就可以了。USACO的解答里甚至只用一个变量就足矣了。具体实施请看代码：

C语言: 高亮代码由发芽网提供

01 #include <stdio.h>
02 #include <stdlib.h>
03 #include <string.h>
04
05 int note[5050];
06
07 int main()
08 {
09     freopen("theme.in","r",stdin);
10     freopen("theme.out","w",stdout);
11     int N,i,d;
12     int len;
13     int maxlen = 0;
14     scanf("%d\n",&N);
15     for(i = 0;i < N;i ++)
16         scanf("%d",&note[i]);
17     for(d = 1;d < N;d ++)
18     {
19         len = 1;
20         for(i = N - d - 2;i >= 0;i --)
21         {
22             if(note[i] - note[i + d] == note[i + 1] - note[i + d + 1])
23                 len ++;
24             else
25                 len = 1;
26             if(len > d)    len = d;   //themes do not overlap, careful here
27             if(maxlen < len)    maxlen = len;   //compare while computing
28         }
29     }
30     if(maxlen >= 5)
31         printf("%d\n",maxlen);
32     else
33         printf("0\n");   //no theme found
34     return 0;
35 }