基数排序——你值得拥有

# 基数排序 ## 引入 　　我们知道“基于比价的排序算法时间复杂度下限是$O(n * log_2n)$”（可以简单地理解为给定序列的排列方式有$n!$种，然后我们每次进行一次比较最多可以排除掉一半的序列，最后复杂度相当于$O(log_2(n!))=O(nlog_2n)$（斯特林公式可求））。 　　但是现实中算法总会有常数，而且快速排序会被卡（C++对sort的优化还可以，例如根据数据特征判断不同时候的方法等等），一般情况下stable_sort采用的是归并排序和插入排序（详见STL源码），其实效率有时还会高一些。 ## 桶排序 　　但是如果希望一个$O(n)$的排序怎么办？我们学习过一个传统的桶排序，这里稍微介绍一下桶排序，时间复杂度为$O(n+max\{val[i]\})$。 > 　　桶排序的原理是按照值把待排序的分类并且存入容器中，然后顺序扫描容器，就可以按照顺序得到所有元素，例如  　　在下一节准备了一份样例代码。 ## 桶排序的多关键字排序 　　桶排序如何做多关键字排序？这对于我们的基数排序算法十分重要。桶排序的多关键字排序主要依赖于反复颠倒排名与位置的关系，也就是这样的一张图：  　　同时建议阅读下面的注释辅助和下面的解释食用。 ```cpp #include <bits/stdc++.h> using namespace std; const int MAXN = 100; const int MAXM = 1000; int a[MAXN], b[MAXN]; int c[MAXM + 1], d[MAXM + 1], f[MAXM + 1]; struct Data{ int x, y; }e[MAXN]; int main(){ srand(233333333); int n; scanf("%d", &n); //生成随机数据 for(int i = 1; i <= n; i++) rand(); for(int i = 1; i <= n; i++) printf("%5d ", a[i] = rand() % MAXM + 1); puts(""); //单关键字排序 for(int i = 1; i <= n; i++) c[a[i]] ++;//桶的基本操作 for(int i = 1; i <= MAXM; i++){ for(int j = 1; j <= c[i]; j++){ printf("%5d ", i); } } puts(""); //单关键字排序输出排名 for(int i = 1; i <= MAXM; i++) c[i] += c[i - 1];//累加排名 for(int i = 1; i <= n; i++) printf("%2d ", c[a[i]] --);//输出的是排名，倒过来就可以了 puts(""); //生成随机数据 for(int i = 1; i <= n; i++) scanf("%d", &e[i].x); // printf("%5d ", e[i].x = rand() % MAXM + 1); puts(""); for(int i = 1; i <= n; i++) scanf("%d", &e[i].y); // printf("%5d ", e[i].y = rand() % MAXM + 1); puts(""); //双关键字排序 memset(c, 0, sizeof(c)); for(int i = 1; i <= n; i++) c[e[i].x] ++; for(int i = 1; i <= MAXM; i++) c[i] += c[i - 1];//c此时记录的是数字为i的数的最高排名 for(int i = 1; i <= n; i++) d[e[i].y] ++; for(int i = 1; i <= MAXM; i++) d[i] += d[i - 1];//d同上面的c for(int i = 1; i <= n; i++) b[d[e[i].y] --] = i;//从优先级比较低的关键字开始，这里记录的是排名为i的e[i].y的位置 for(int i = n; i >= 1; i--)//我们如何处理第一关键字相等的情况？如果第一关键字相等的话，那么第二关键字在前面的会被后访问到（循环倒序），同时第一次访问同一个e[]的时候的时候c数组--，然后因此后第二个访问的位置由一个排名更前的c（c会更小）去填充 a[c[e[b[i]].x] --] = b[i];//a这里记录的是排名为i的e[i].x的位置 for(int i = 1; i <= n; i++) printf("%2d ", a[i]);//如果我们进行一次c[a[i]] = i的这样的变换，就可以变换排名与位置的关系 puts(""); for(int i = 1; i <= n; i++) printf("%5d ", e[a[i]].x); puts(""); for(int i = 1; i <= n; i++) printf("%5d ", e[a[i]].y); puts(""); return 0; } ``` 　　为什么`d[e[i].y]] --`可以取得排名对应的呢？首先考虑到我前面统计的d长什么样子。首先d作为桶有值的地方都++，然后再统计了一发前缀和，于是**得到的就是在没有重复的情况下得到的排名**，然后每次遇到一次`e[i].y`的时候`d--`，然后就把同一个数分配了几个不同的排名，例如我第一次遇到9的时候排名是8，第二次遇到的时候就是7了。但是注意我们分配排名的时候是递减的 　　然后第二次`a[c[e[b[i]].x] --] = b[i]`的时候我们的循环**倒序**枚举的是排名，首先最里层将取出排名为i的第二关键字的位置，然后向外一层取出该位置的第一关键字的值，然后再向外一层取出第一关键字的值的排名，同时注意到第一关键字相同但是第二关键字靠后的将会在靠前的次序取出，这个时候**分配的第一关键字靠后**，于是完成了排序。 　　对于多层排序的分析也就像上面一样一层层解析数组即可。 ## 基数排序 　　我们发现，如果我有一个非常大的元素需要排序的话那么这个桶是存不下的，同时最大的元素大了之后扫描桶的时间复杂度非常高。 　　所以如何优化这个东西呢？于是就有基数排序。 　　标准的基数排序每次从小到大排序十进制下的一个位，**并且只关心这一个位**，例如下图  　　为什么要从低位到高位排：首先是为了保证相同高位的数字可以能够按照低位排序，而且如果高位优先，那么高位排好的序在低位又会被打乱。 ## 优化的基数排序 　　注意到上面操作用的时间是$(n + log_{10}(Max\{val[i]\})) * (log_{10}(Max\{val[i]\}))$，简写为$O(nlog_rm)$，其中$n$为元素个数,$r$是排序所用的基数，$m$为所有元素最大值的大小。 　　但是一般的$log_rm$接近10，而且还涉及到取模等操作，速度略慢，常数较大。 　　所以我们考虑一下以二的幂次为基数的排序，首先一个比较好的优点应该是它的几乎所有操作都可以围绕二进制操作，常数小，然后再一个就是如果以$2^8=256=0×ff+1$为一个基数的话,那么int范围内的$log_rm$只有4，然后255的数组还可能可以写入内存而不是硬盘（也有可能写入缓存）。 　　于是我们发现以二的幂次为基数的二进制基数排序是我们所需要的，而且效果也比较优秀，至少目前[Luogu1177快速排序](https://www.luogu.org/problemnew/show/P1177)的某个匿名者用的是这个优化，效果便是目前的16ms的排序（当然我们可以确信它可以做得更好）。 　　这张图描述了一次这样的排序。  ## 一个扩展的功能 　　同时我们有时候还会有一个对于结构体进行多关键字排序的需求，我们注意到C++内定的数据类型至少为1个字节，一般用的整数为4个字节，然后我们发现对于一个32为整数$k$，我们对k进行32为的比较可以认为是先进行16位高位的比较然后再进行16位的低位的比较，在高位相等的时候比较低位，这不就是双关键字排序吗？ 　　于是结构体的双关键字排序可以这样做：先利用指针类型转换结构体成为32/64位整数，然后再适当调整基数直接排序这些整数，然后就可以得到这些排序后的结构体。 　　结构体的内存分布如下：  　　对于这个结构体，使用64位排序，基数$2^{16}$（$2^8$也行） ## 实现方法 　　这里直接以template作为参数进行32位排序，这个参数接受不高于32位的类型。a为待排序数组，b为临时数组 　　现在看看代码如何写。 1. 首先开四个数组作为4个关键字的桶（基数为$2^8$） 2. 把元素用二进制运算拆开，放入桶内 3. 桶子进行桶排序，利用上一次排序的结果交替利用数组进行下一次排序 　　注意结构体排序的待排序关键字应该定义在结构体的前面（和内存有关），不建议使用union。同时在以前所有r数组开的都是0xff的大小,但是后来考虑到实际使用的时候不同的机子或不同的时间下有不同的内存情况,所以后来把数组扩大了1(测试过不会导致缓存的问题). ```cpp template <typename T> inline void Radix_Sort(register const int n, T *a, T *b){//模板匹配参数 int r1[0x100], r2[0x100], r3[0x100], r4[0x100]; memset(r1, 0, sizeof(r1)), memset(r2, 0, sizeof(r2)); memset(r3, 0, sizeof(r3)), memset(r4, 0, sizeof(r4)); register int i, tmp_int; register T * j, *tar;//注意C++11 下已经基本弃用register，所以不建议使用C++11去尝试卡常（C++下优化明显） for(j = a + 1, tar = a + 1 + n; j != tar; ++j){//这一部分用的是指针加速判断（O2下优化比较大） tmp_int = * (int *) j;//用指针转换类型 ++r1[tmp_int & 0xff];//二进制取低8位 ++r2[(tmp_int >> 8) & 0xff];//二进制取16位 ++r3[(tmp_int >> 16) & 0xff];//24位 ++r4[tmp_int >> 24];//32位 } for(i = 1; i <= 0xff; ++i){//桶累加统计 r1[i] += r1[i - 1]; r2[i] += r2[i - 1]; r3[i] += r3[i - 1]; r4[i] += r4[i - 1]; } for(j = a + n; j != a; --j){//排名统计 tmp_int = * (int *) j; b[r1[tmp_int & 0xff]--] = *j;//二进制避免了取模 } for(j = b + n; j != b; --j){//看上去细节不少，但是其实发现这几行代码都可以复制 tmp_int = * (int *) j; a[r2[(tmp_int >> 8) & 0xff]--] = *j; } for(j = a + n; j != a; --j){//参考桶排序的代码去理解 tmp_int = * (int *) j; b[r3[(tmp_int >> 16) & 0xff]--] = *j; } for(j = b + n; j != b; --j){ tmp_int = * (int *) j; a[r4[tmp_int >> 24]--] = *j; } } ``` ## 此文章还带有更新 ### 更新1 ~~好消息好消息，浮点数可以基数排序了~~ ~~好消息好消息，任意大小结构体可以基数排序了~~ 我们需要感谢[IEEE 754](https://en.wikipedia.org/wiki/IEEE_754)的标准，它制定的浮点标准还带有一个特性：如果不考虑符号位（浮点数是按照类似于原码使用一个bit表示符号位s，s = 1表示负数），那么同符号浮点数的比较可以等价于无符号整数的比较，方法就是先做一次基数排序，然后再调整正负数部分顺序即可 同时还提供一种不止64bit的待排序位数的结构体排序，代码在下面，会被long double的排序调用到，但是特别注意这种方法只能在小端法电脑上用，大端法需要把k的循环反过来 代码如下 ``` cpp #include <cstdio> #include <cstring> #include <algorithm> using namespace std; const int MAXA = 1000001; template <typename T> inline void Radix_Sort_Bit(register const int n, T *a, T *b){ size_t size_of_type = sizeof(T); size_t num_of_buc = size_of_type >> 1; unsigned ** r = new unsigned *[num_of_buc]; register int i, k;// 可以直接开成r[num_of_buc][0x10000] for(i = 0; i < num_of_buc; i++) r[i] = new unsigned [0x10000], memset(r[i], 0, 0x10000 * sizeof(unsigned)); register unsigned short tmp_us; register T * j, *tar; for(k = 0; k < num_of_buc; ++k){// 填充桶 for(j = a + 1, tar = a + 1 + n; j != tar; ++ j){ tmp_us = * (((unsigned short *)j) + k); ++ r[k][tmp_us];// 直接强制转换指针然后计算地址 } } for(k = 0; k < num_of_buc; ++k)// 累加前缀和 for(i = 1; i <= 0xffff; ++i) r[k][i] += r[k][i - 1]; for(k = 0; k < num_of_buc; k += 0x2){// 反复“倒桶” i = k; for(j = a + n; j != a; --j){ tmp_us = * (((unsigned short *)j) + i); b[r[i][tmp_us]--] = *j; } i |= 1; if(i == num_of_buc) break; for(j = b + n; j != b; --j){ tmp_us = * (((unsigned short *)j) + i); a[r[i][tmp_us]--] = *j; } } for(int i = 0; i < num_of_buc; i++) delete[] r[i]; delete [] r; } inline void Radix_Sort_LDouble(register const int n, long double *a, long double *b){ Radix_Sort_Bit(n, a, b); reverse(a + 1, a + 1 + n);// 转换负数顺序 reverse(upper_bound(a + 1, a + 1 + n, (long double)(-0.0)), a + 1 + n); } int n; long double a[MAXA], b[MAXA]; int main(){ freopen("fl.in", "r", stdin); freopen("fl.out", "w", stdout); scanf("%d", &n); for(int i = 1; i <= n; i++) scanf("%Lf", a + i); Radix_Sort_LDouble(n, a, b); for(int i = 1; i <= n; i++) printf("%.18Lf\n", a[i]); fclose(stdin), fclose(stdout); return 0; } ``` 实测会比sort快，但是上面写的是兼容模式下的，更好的实现应该要用int128或者是做循环展开或者对于较小的数据类型把数据开到栈内，建议使用的时候自己按照数据类型调整一下 ## 使用实例 　　[Luogu P2672 推销员](https://www.luogu.org/problemnew/show/P2672) 　　这是一道贪心题，思路比较微妙： > $\ \ \ \ $注意到$a$对于答案的贡献是$\Sigma{a}$，但是$s$对答案的贡献是$2 * Max\{s\}$ > $\ \ \ \ $分析对于推销$k$个人答案的组成：要么答案就是由$a$值最大的$k$个人构成，要么答案是由$a$值最大的$k-1$个人和一个$s$非常大的人构成 　　这里对于排序的需求是对于结构体储存的a,s基于a的单关键字排序。 　　然后这里使用基数排序优化常数，目前Rank1。 ```cpp #include <cstdio> #include <fcntl.h> #include <unistd.h> #include <sys/mman.h> //User's Lib #include <cstring> // #include <algorithm> using namespace std; char *pc; char outp[1111111], *op = outp; inline void Main_Init(){ static bool inited = false; if(inited){ fwrite(outp, 1, op - outp - 1, stdout); fclose(stdin), fclose(stdout); } else { #ifndef ONLINE_JUDGE freopen("b.in", "r", stdin); #endif pc = (char *) mmap(NULL, lseek(0, 0, SEEK_END), PROT_READ, MAP_PRIVATE, 0, 0); inited = true; } } inline int read(){ int num = 0; char c; while((c = *pc ++) < 48); while(num = num * 10 + c - 48, (c = *pc ++) >= 48); return num; } inline void Write(const int &tar){ if(!tar) return ; Write(tar / 10); *op ++ = (tar % 10) ^ 48; } namespace LKF{ template <typename T> extern inline T abs(T tar){ return tar < 0 ? -tar : tar; } template <typename T> extern inline void swap(T &a, T &b){ T t = a; a = b; b = t; } template <typename T> extern inline void upmax(T &x, const T &y){ if(x < y) x = y; } template <typename T> extern inline void upmin(T &x, const T &y){ if(x > y) x = y; } template <typename T> extern inline T max(T a, T b){ return a > b ? a : b; } template <typename T> extern inline T min(T a, T b){ return a < b ? a : b; } } //Do not copy from other files //Source Code const int MAXN = 100011; struct Data{//第一关键字a在前面，占32bit int a, s; }data[MAXN], tmp_data[MAXN]; template <typename T> inline void Radix_Sort(register const int n, T *a, T *b){ int r1[0x100], r2[0x100], r3[0x100], r4[0x100]; memset(r1, 0, sizeof(r1)), memset(r2, 0, sizeof(r2)); memset(r3, 0, sizeof(r3)), memset(r4, 0, sizeof(r4)); register int i, tmp_int; register T * j, *tar;//让编译器辅助指针优化 for(j = a + 1, tar = a + 1 + n; j != tar; ++j){ tmp_int = * (int *) j; ++r1[tmp_int & 0xff]; ++r2[(tmp_int >> 8) & 0xff]; ++r3[(tmp_int >> 16) & 0xff]; ++r4[tmp_int >> 24]; } for(i = 1; i <= 0xff; ++i){ r1[i] += r1[i - 1]; r2[i] += r2[i - 1]; r3[i] += r3[i - 1]; r4[i] += r4[i - 1]; } for(j = a + n; j != a; --j){ tmp_int = * (int *) j; b[r1[tmp_int & 0xff]--] = *j; } for(j = b + n; j != b; --j){ tmp_int = * (int *) j; a[r2[(tmp_int >> 8) & 0xff]--] = *j; } for(j = a + n; j != a; --j){ tmp_int = * (int *) j; b[r3[(tmp_int >> 16) & 0xff]--] = *j; } for(j = b + n; j != b; --j){ tmp_int = * (int *) j; a[r4[tmp_int >> 24]--] = *j; } } int val[MAXN]; int main(){ Main_Init();//如果怀疑你的结构体的内存排布不对，可以用sizeof(Struct)查看占用字节数 int n = read(); for(register int i = 1; i <= n; ++ i) data[i].s = read(); for(register int i = 1; i <= n; ++ i) data[i].a = read(); Radix_Sort(n, data, tmp_data); for(register int i = 1, j = n; i < j; ++ i, -- j) LKF::swap(data[i], data[j]); for(register int i = n; i >= 1; i--){ val[i] = LKF::max(val[i + 1], (data[i].s << 1) + data[i].a);//下面这个奇怪的函数是优化缓存，也就是提前把一个函数扔到高速缓存内以供CPU使用。第一个参数是目标需要写入缓存的变量的指针。第二个参数表示读(0)还是写(1)，第三个参数表示时间局部性，时间局部低的变量用过一次之后很长一段时间内不会使用，也就是在缓存内时间小，该参数从0到3表示时间局部性越来越高 __builtin_prefetch(val + i - 1, 1, 0); } int Sum1 = 0, Sum2 = 0, Max = 0; for(register int i = 1; i <= n; ++ i){ Sum1 += data[i].a; LKF::upmax(Max, data[i].s); Write(LKF::max(Sum1 + (Max << 1), Sum2 + val[i])), *op ++ = '\n';//两种答案取最大 Sum2 += data[i].a; __builtin_prefetch(data + i + 1, 0, 0); } Main_Init(); return 0; } ``` ## 课后思考 　　我在上一份代码中间留下了一个瑕疵：我首先对data数组进行了从小到大的排序，然后再通过swap调换顺序。那么如何修改基数排序使其资瓷从大到小的排序呢？请思考。 ## 写在最后 　　其实学习基数排序会比较实用，至少理解了这个之后就可以对后缀数组的倍增求法的理解有大大的长进。 　　然后绘图工具：GraphViz和Sequence Diagrams。 　　参考不少博客（在此不一一列出）和《算法导论》对基于比较的排序算法的下限的分析。