算法

二分搜索虽然基本思想简单，但其细节却令人意外的抓狂（Although the basic idea of binary search is comparatively straightforward, the details can be surprisingly tricky）。这里我们来分析一下常见写法的坑。 两类不同的二分 二分有两类，一是找值，即判断某个值存不存在，二是找边界，前者比后者简单很多，以下是找值的典型写法（来自中文wiki百科） int binary_search(const int arr[], int start, int end, int key) { int ret = -1; // 未搜索到数据返回-1下标 int mid; while (start <= end) { mid = start + (end - start) / 2; //直接平均可能會溢位，所以用此算法 if (arr[mid] < key) start = mid + 1; else if (arr[mid] > key) end = mid - 1; else { // 最後檢測相等是因為多數搜尋狀況不是大於要不就小於 ret = mid; break; } } return ret; // 单一出口 } 对于找边界，有四种情况，如以下例子，查找数值5的四种边界 1 1 1 2 2 5 5 5 5 7 9 ^ 1 小于的最右元素 ^ 2 大于等于的最左元素 ^ 3 小于等于的最右元素 ^ 4 大于的最左元素 既然有四种边界，那就有四种写法。

在上一篇介绍了基础优化算法后，本篇介绍更复杂的内容。本篇的三大内容：FFT，NTT，分治除法。中文资料中我尚未发现有博客文章在介绍分治除法的，所以我就来写第一个介绍吧。 FFT (Schönhage–Strassen algorithm) FFT就是快速傅里叶变换的缩写，FFT这里不重点介绍，参见oiwiki中对FFT的介绍 用比较显浅且简要的方式来介绍的话，FFT可以在$O(nlogn)$时间里（注：严格来说应该是$O(nlognloglogn)$，论文链接在这，在274页7.8处开始证明）求出两个多项式的乘积，而大整数乘法就可以视为两个整数系数的多项式的积再处理进位即可。那这样就产生一个问题，如果这个大整数采用的基数base，且长度为n，那么多项式的积的最大值可能会达到$n(base-1)^2$，因为FFT采用浮点运算，那么这个最大值必须在double的精度内，所以如果采用万进制，那么很容易导致精度不足进而结果错误，另外，还有求单位复根时的误差也会让它提早出现结果错误。实际应用的时候，一般采用10进制来增加长度的上限。所以如果你的大整数采用压位表示，那么在运用FFT之前，要先做拆位拆成10进制。总之，因为浮点精度的问题，不建议在大整数实现里使用FFT，建议采用下文要介绍的NTT来代替。所以这里对FFT的优化不做介绍，重心放在NTT上。 NTT NTT就是快速数论变换，它和FFT很像，但它在整数域上做变换，具体数学介绍参见oiwiki中对NTT的介绍。由于没有浮点数参与，于是没有精度上的问题。

在上一篇介绍了基础算法后，本篇介绍算法级别的优化。本篇的三大内容：乘法优化，任意进制输入，任意进制输出。 乘法优化 这里要介绍的，是Karatsuba发现的分治法，我们假设要相乘的两个数，都有2n位，那么这两个数就可以分别表示为$a_1base^n+a_2, b_1base^n+b_2$，其中，$a_1,a_2,b_1,b_2$是n位的大整数，那么，它们的积就是 $$\begin{align} & (a_1base^n+a_2) \times (b_1base^n+b_2) \\ & = a_1b_1base^{2n} + (a_1b_2+a_2b_1)base^n + a_2b_2 \\ & = a_1b_1base^{2n} + ((a_1+a_2)(b_1+b_2)-a_1b_1-a_2b_2)base^n + a_2b_2 \end{align}$$ 如果这样不够明显的话，我们用$c_1$代替$a_1b_1$，用$c_3$代替$a_2b_2$，得到 $c_1base^{2n} + ((a_1+a_2)(b_1+b_2)-c_1-c_3)base^n + c_3$ 于是，这里一共有3次乘法，比起原来的4次暴力乘法减少了1次。而里面的乘法又可以进行递归优化，时间复杂度从$O(n^2)$下降到$O(n^{log_23})$约$O(n^{1.585})$

最近在编写大整数库的过程中，踩到不少的坑，于是把一些有用的细节准备写成文章做整理。如果你只是想直接查找并使用一个大整数库，那直接上GMP即可，如果是用在比赛，那直接用我的MiniBigInteger项目，如果你是想学习个中细节，那你可以坐下来细品。 所谓大整数，又叫高精度运算，就是运算对象是上千位甚至到百万位，总之远远超过内置数据类型的表示范围，这类数字都叫大整数。而C/C++的标准库里目前并没有大整数库，于是这个轮子被反复制造了无数个，不过在github上比较有质量的轮子并没有很多。本文除了介绍基础实现，主要还是介绍优化方法。 大整数的表示 大整数的表示方法最常见的有4种： 直接使用string 使用定长数组(仅适用于竞赛) 使用链表 使用变长线性表例如vector 直接用string的方式适合初学者，输入输出直观，但缺点也非常明显，因为计算时需要在字符与数值之间来回转换，浪费太多不必要的时间，效率会非常差。不过如果你是初学者，先用string表示法来写未尝不是个好主意。但有个细节就是，如果想要效率高，最好把string前后倒置调整为低位在前再做运算，这样速度和实现难度都会低一些。 至于使用链表，好处是变长容易，变短也不难，但性能比用string的更差还更难写，这里就不谈了，以下介绍使用数组的表示法 为了在数组里表示一个大整数，如果我们采用10进制，表示123456789，那很简单，例如这样：

这次主要介绍伪随机数生成算法，顺便介绍一个在2018-2019年的伪随机数研究成果，就是 xoshiro/xoroshiro 随机数生成算法。 历史 在较早的时候，甚至到现在，伪随机数的生成元老级别算法“线性同余伪随机数生成算法”可谓无处不在，像现在的C/C++的rand函数就是使用线性同余实现的伪随机数生成。所谓的线性同余法，就是这个迭代方程 $S_n = (aS_{n-1} + c)\mod m$，其中，$S_0$ 称为这个随机序列的种子，a,c,m是三个常数，不过这三个数不能随意选，m的大小决定随机数的周期，最大周期等于m，为了便于在计算机里实现，通常m选取$2^{32}$或$2^{64}$。在m已经确定为这两的时候，为了让周期尽可能大，常数a,c还至少要满足以下条件： 若 c 非 0，那么 c 与 m 互质；若 c 为 0，那么 a 与 m 互质 m 所有的素因子均能整除 a-1 若 m 是4的倍数，那么 a-1 也是 4 的倍数 a 和 c 都是正整数且小于 m 一些典型的常数取值可以参见wiki上的线性同余条目。 更高的需求 我们之所以使用线性同余，就是因为它实现简单，在对随机数质量要求较低的时候，例如用来作为treap的随机数，那么线性同余完全够用，但建议不要使用rand，因为在windows下不少编译器的最大值太小了，导致效果下降，自己写一个用参数a=69069,c=1,m=2^32比rand好，我在那篇关于treap的文章就是用了这组参数。线性同余法最大的缺陷是低位随机性特别差，如果使用类似next() % k的方式来获得区间在$[0,k-1]$的随机数，那么当线性同余迭代方程的m是2的幂且k也是2的幂的时候，灾难就发生了，特别地当k是2的时候，你将得到一个0101的循环序列。为了避免这种情况，通常会取线性同余结果的高位，而且低位去掉得越多，%2的周期就越长。例如结果是64位，取高32位作为最终结果，那么%2的周期就是 $2^{33}$ ，但这样会导致有效位减少，而且问题也没有根本地解决。另一种解决办法是选取一个素数作为m，例如2147483647正是一个素数，但如此一来，线性同余的计算速度就会慢不少，周期也没有前一个的长。两种基本实现如下：

KMP已经介绍过了，这次主要介绍的是Manacher 最长回文子串 Manacher算法要解决的问题就是求最长回文子串，用到的思维和扩展KMP实在是像，不过理解起来比扩展KMP简单。 先定义数组v，v[i]表示以第i个字符为中心，到回文串一端的距离，我们以字符串"acabacab"为例，如下表（index是下标） string a c a b a c a b \0 index 0 1 2 3 4 5 6 7 8 m ^ i ^ e ^ v 1 2 1 4 i是当前要计算的指针，m是上次计算的指针，e是下一个要比较的位置的指针 然后++i，注意这时候，由于以字符b两边是对称的，所以在求v[4]的值的时候，可以先查v[2]的值，是1，所以v[4]也是1 string a c a b a c a b \0 index 0 1 2 3 4 5 6 7 8 m ^ i ^ e ^ v 1 2 1 4 1 再继续++i，同样的，在求v[5]的值的时候，可以先查v[1]的值，是2，但是，这个长度达到了e指针的位置，即i+2>=e，这时候就更新指针m，并扩展e的位置，即比较str[7]与str[3]，找到以下标5为中心的回文串边界。

KMP之所以在竞赛中常见，并不是因为它用来匹配字符串，而是用它的next数组，为了介绍它，我们先讲讲最长公共前缀 最长公共前缀 我们拿字符串ababcab作为例子 string a b a b c a b len 0 0 1 2 0 1 2 这里所表达的是，例如取第3、4个字符”ab”，这个子串与前缀完全匹配，且它的长度是2，所以就记录2，而第3、4、5个字符”abc”与前缀不能完全匹配，就记作0，含义就这么简单，而且你会发现，计算b的时候，可以根据它所匹配的字符的偏移来，b如果是匹配的，就找到匹配的那个字符是数组中的第几个，它是第二个，所以填2进去。我们再来看更复杂的例子 string a b a c a b a b len 0 0 1 0 1 2 3 2 最后那个字符不匹配的时候，1是怎么计算出来的呢，直接重新计算当然也可以，但就出现重复计算了。我们考虑一下匹配过程，在前面的字符a的时候，前后各一个指针，像这样 string a b a c a b a b len 0 0 1 0 1 2 3 ? pointer ^ ^ 然后两个a匹配，arr[6] = pointer1 - arr 得到3，然后两指针一起移动 string a b a c a b a b len 0 0 1 0 1 2 3 ? pointer * ^ ^ 这时候，不匹配，那么前一个指针上一次指向的是arr[2]的位置，即图上*的地方，值是1，这个值如果是p，那就移动到arr[p]的地方，所以就移动到arr[1]的地方，本质上就是找到前一个匹配此后缀的位置，即 string a b a c a b a b len 0 0 1 0 1 2 3 2 pointer ^ ^ 然后再尝试匹配，这次匹配上了，然后前一指针指向第二个元素，所以赋值2

在上一篇我们介绍了快排的各种优化，最后得到了一个超越VS的STL std::sort实现的版本，这回我们来讲讲怎么折磨gcc的std::sort。 GCC的实现 这个我直接给个github来源，链接里是gcc8分支的实现源码，我已经通过链接标记出第1896行，那里就是__unguarded_partition函数的实现，就是快排的划分函数，而在后面几个函数就是快排的本体了。为了避免它代码更新导致位置变化，我把相关代码复制过来。

在上一篇我们介绍了四种不同的划分算法，现在我们来针对Hoare partition scheme来讲解一些优化和注意的问题。 最坏时间复杂度的优化 在前一篇的示例代码里面，只是最简单地选择的最开头或最后面的元素作为划分，这对于乱序的数据还好，对于有序的数据这么做，时间复杂度就直接变成$O(n^2)$了，那么怎么解决？第一个要解决的反而不是划分元素的选择上，而是改成intro sort，记录递归深度或类似的办法，到达限制条件的时候改而使用堆排序，这属于混合排序，先让排序的最坏时间复杂度降下来是第一要务。所以可以改写出以下代码：

本篇补充smooth sort和weakheap sort，不打算做太多介绍，因为自己太菜不会讲，只作为笔记来记录一下实现代码。smooth sort在接近排序完成的情况下的动态图