【学习笔记】我们仍未理解那天所学习的数论算法

前言

对于部分基础的数论知识点（如整除、同余等），笔者不会再进行讲解。因此，笔者希望有扎实的数论基础的读者来阅读这篇文章。

中国剩余定理

对于线性同余方程组 $x \equiv a_i \pmod{p_i}$ ，中国剩余定理能够求解其最小非负整数解。

这个算法要分为两种情况讨论，而两种情况的算法实际上存在本质区别。

模数两两互质

我们有如下基于构造思想的算法：

设 $\displaystyle \text{M} = \prod_{i = 1}^mp_i$ ， $\text{M}_i = \dfrac{\text{M}}{p_i}$ ， $\text{M}_i^{-1}$ 为 $\text{M}_i$ 在模 $p_i$ 意义下的逆元，那么整个方程组的最小非负整数解为 $\displaystyle (\sum_{i = 1}^m a_i\cdot \text{M}_i\cdot \text{M}_i^{-1}) \bmod \text{M}$ 。

证明：

考虑这个和式的第 $i$ 项 $a_i \cdot \text{M}_i \cdot \text{M}_i^{-1}$ ，将其对 $p_i$ 取模即得到第 $i$ 个同余方程的特解 $a_i$ ，而将其对其它的 $p_j~(j \ne i)$ 取模均得到 $0$ ，所以这个项既能够满足第 $i$ 个同余方程，又不会对其它的同余方程造成影响。

这样，我们证明了这个和式是原方程的一个特解。

又因为方程组在 $[0, \text{M})$ 内仅存在一个解，因此将这个和式对 $\text{M}$ 取模即得到最小非负整数解。

值得一提的是，我们熟知的拉格朗日差值实际上是基于这一部分的算法的。

模数不两两互质

我们是否能够考虑套用上述算法的框架呢？

不行，因为 $\text{M}_i^{-1}$ 甚至有可能根本不存在。因此，仅对上述算法进行微调是不可行的。

我们考虑一个很简单的思想：合并同余方程。

对于前 $i$ 个同余方程，我们可以通过将它们合并得到一个形如 $x \equiv b_i \pmod{\text{M}_i}$ 的方程，这里的 $\text{M}_i$ 表示前 $i$ 个模数的 $\operatorname{lcm}$ 。

接下来我们考虑其和第 $i + 1$ 个方程的合并。

首先，我们将满足前 $i$ 个方程的 $x$ 写成 $b_i + \text{M}_i \cdot k~(k \in \Z)$ 的形式，那么我们仅需求出 $k$ 的一个特解便可得到前 $i + 1$ 个方程的特解。然后，我们将其带入第 $i + 1$ 个方程得到 $b_i +\text{M}_i \cdot k \equiv a_i \pmod{p_i}$ ，再将其转化为二元一次不定方程的形式，最后用 $\operatorname{exgcd}$ 求出 $k$ 的一个特解即可。在这个过程中，我们还可以判断方程组是否存在无解的情况。

实际上，若我们对于每一个方程，都将 $x$ 带上一个 $c_i$ 的系数，那么无论是上述的哪一种情况，对应的算法经过调整仍然是可行的。

读者不妨自行思考，在 $x$ 带上了系数的情况下，两种情况对应的算法分别应该如何调整。

Lucas 定理

我们经常需要计算 $\displaystyle \binom{n}{m} \bmod p$ 的值，而若 $n$ 和 $m$ 较大，传统的算法是无法支持的。

是否有其它的算法呢？

模数为质数

$\mathcal{Lucas}$ 定理是这样一个公式：

$\binom{n}{m} \equiv \binom{\displaystyle \lfloor \frac{n}{p} \rfloor}{\displaystyle \lfloor \frac{m}{p} \rfloor} \cdot \binom{n\bmod p}{m\bmod p} \pmod{p}$

根据这个公式，我们可以预处理 $n, m \in [0, p)$ 的答案，然后便可做到单次时间复杂度 $\mathcal{O}(\log_pn)$ 地计算组合数。若 $p$ 的范围较小，那么这是一个复杂度十分优秀的算法。

$\mathcal{Lucas}$ 定理还有一个十分优美的推论：

$\displaystyle \binom{n}{m} \bmod 2$ 的值为 $1$ ，当且仅当 $n~\texttt{and}~m = m$ ，也即在二进制下 $m$ 是 $n$ 的子集，其中 $\texttt{and}$ 表示二进制下的按位与运算。

因为当 $p = 2$ 时， $\mathcal{Lucas}$ 定理会将 $\displaystyle \binom{n}{m}$ 拆分为若干个 $\displaystyle \binom{0}{0},\binom{0}{1},\binom{1}{0},\binom{1}{1}$ ，而其中只有 $\displaystyle \binom{0}{1}$ 等于 $0$ ，这正好对应了 $n~\texttt{and}~m \neq m$ 的情况。

这个推论揭示了组合数的奇偶性与 $n$ 和 $m$ 的二进制之间的关系。对于和组合数的奇偶性相关的问题，将其转化为和 $n$ 和 $m$ 的二进制相关的问题，不失为一种好的办法。

模数不为质数

当 $p$ 不为质数时，上述定理是不成立的，因为其证明用到了和质数相关的性质。

对于模数为合数的求值问题，将其模数进行质因数分解，最后再用中国剩余定理合并答案，不失为一种好的办法。

我们按照上述方法操作，将问题转化为如何快速求 $\displaystyle \binom{n}{m} \bmod p^q$ 。

将 $\displaystyle \binom{n}{m}$ 拆开得到 $\displaystyle \frac{n!}{m! \cdot (n - m)!}$ ，然后把三个部分包含的 $p$ 全部提取出来，原式将会转化为 $\displaystyle \frac{\displaystyle \frac{n!}{p^{r_1}}}{\displaystyle \frac{m!}{p^{r_2}} \cdot \frac{(n - m)!}{p^{r_3}}} \cdot p^{r_1 - r_2 - r_3}\bmod p^q$ ，此时分式的分母与模数互质，那么问题的关键就在于如何快速求 $n! \bmod p^q$ 。

同 $\mathcal{OI-Wiki}$ ，这里用 $22! \bmod 3^2$ 来举例。

先把 $22!$ 拆开：

$22! = 1\times2\times3\times4\times5\times6\times7\times8\times9\times10\times11\times12\times13\times14\times15\times16\times17\times18\times19\times20\times21\times22$

考虑把 $3$ 的倍数提出来，化简一下就有：

$22! = 3^7\times(1\times2\times3\times4\times5\times6\times7) \times(1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16\times17) \times (19\times20\times22)$

第一个部分是 $p$ 的 $\displaystyle \lfloor\frac{n}{p}\rfloor$ 次幂，我们需要记录其指数以计算 $r$ 。

第二个部分是 $\lfloor\dfrac{n}{p}\rfloor!$ ，我们可以递归求解。

第三个部分是 $1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16\times17$ ，仔细观察可以发现：

$1\times2\times4\times5\times7\times8 \equiv10\times11\times13\times14\times16\times17\pmod{3^2}$

形式化地，有：

$\prod_{i = 1,\gcd(i, p) = 1}^{p^q} i \equiv \prod_{i = 1,\gcd(i, p) = 1}^{p^q} (i + p^q \cdot s) \pmod{p^q}$

（其中 $s \in \Z$ 。）

这个循环节一共循环了 $\lfloor\dfrac{n}{p^q}\rfloor$ 次，我们直接将其求出来，然后再对其求幂即可。

最后是 $19 \times 20\times 22$ ，形式化的表示即为：

$\prod_{i = 1,\gcd(i, p) = 1}^{n\bmod p^q} i$

我们在求解第三个部分的答案时顺便计算其值即可。

实际上，如果询问次数较多，我们还可以对第三个部分和第四个部分进行预处理，这样单次时间复杂度就仅和 $p^q$ 的个数、第一个部分和第二个部分有关了。

阶与原根

在很多和数论相关的证明中，我们都能看到原根的影子。

这是因为原根所具有的优美性质：

在模 $p$ 意义下，任意一个与 $p$ 互质的数都能够被表示为 $g^i$ 的形式（其中 $g$ 为模 $p$ 的原根， $a\in[0, \varphi(p))$ ），换言之， $g^i$ 遍历整个模 $p$ 的缩系。

接下来我们一起来探讨原根所具有的优美性质。

（注意：下文将出现较多的定理，笔者不会对它们进行证明，对证明感兴趣的读者可以参考 $\mathcal{OI-Wiki}$ 。）

阶

在了解原根之前，我们首先需要了解阶。

设 $a，p$ 互质且 $p > 1$ ，则使得同余式 $a^n \equiv 1\pmod{p}$ 成立的最小的正整数 $n$ 被称为 $a$ 在模 $p$ 意义下的阶，记为 $\delta_p(a)$ 。

我们不妨先来探讨阶的性质。

首先是两个较为显然的性质：

$a, a^2, a^3, \cdots a^{\delta_p(a)}$ 在模 $p$ 意义下两两不同余。

若存在 $n$ 使得 $a^n \equiv 1\pmod{p}$ ，那么有 $\delta_p(a) \mid n$ 。

两条性质均可根据阶的最小性反证。

接下来是两个与四则运算相关的性质：

设 $p \in \N_+$ 和 $a, b \in \Z$ ，且 $\gcd(a, p) = \gcd(b, p) = 1$ ，那么 $\delta_p(a \cdot b) = \delta_p(a) \cdot \delta_p(b)$ 和 $\gcd(\delta_p(a), \delta_p(b)) = 1$ 互为充要条件。

设 $p \in \N_+$ 和 $a, i \in \N$ ，且 $\gcd(a, p) = 1$ ，则有 $\delta_p(a^i) = \dfrac{\delta_p(a)}{\gcd(\delta_p(a), i)}$ 。

两条性质把阶与乘除法运算联系了起来，方便了我们对阶的转化，对定理的证明等都有很大帮助。

原根

若存在 $g \perp p$ 使得 $\delta_p(g) = \varphi(p)$ ，则称 $g$ 为模 $p$ 的原根。

这个定义能够直接证明上文所提到的 “ $g^i$ 遍历整个模 $p$ 的缩系” 这个性质。

我们在探究原根前，首先需要知道一个数是否存在原根。

原根存在定理：正整数 $n$ 存在原根，当且仅当 $n = 1, 2,4,p^q, p^q\cdot 2$ ，其中 $p$ 为奇质数， $q$ 为正整数。

下面我们来挖掘原根的一些性质。

若 $p$ 存在原根，其个数为 $\varphi(\varphi(p))$ 。

原根判定定理：设 $p$ 存在原根，且有 $\gcd(g, p) = 1$ ，那么 $g$ 为模 $p$ 的原根的充要条件是，对于 $\varphi(p)$ 的每个质因数 $d$ ，都不满足 $g^{\frac{\varphi(p)}{d}} \equiv 1 \pmod{p}$ 。

其中第二条性质尤为重要，它使我们能够在对 $\varphi(p)$ 进行质因数分解后，用单次 $\mathcal{O}(\log_2\varphi(p))$ 的时间复杂度判定原根。

那么，对于一个存在原根的数 $p$ ，我们该如何求出其所有原根呢？

我们先考虑枚举出 $p$ 的最小原根 $g$ 。这一步的复杂度基于我国数学家王元于 $1959$ 年的证明：若 $p$ 存在原根，那么最小原根是 $\mathcal{O}(p^{\frac{1}{4}})$ 级别的。所以暴力枚举的时间复杂度在可承受范围内。

枚举出最小原根之后，我们可以枚举所有与 $\varphi(p)$ 互质的 $i \in [1, \varphi(p)]$ ，那么 $g^i \bmod p$ 均为模 $p$ 的原根。该算法的原理可以参考阶的第四条性质，同时该算法的成立也直接证明了 “若 $p$ 存在原根，其个数为 $\varphi(\varphi(p))$ ” 这一条行之。

可以发现，原根的性质大多都能为优化和其相关的算法提供帮助，这便是其的优美之处之一。

熟练运用好原根的性质，对数论题（尤其是数论题的证明）有很大帮助。

离散对数

离散对数实际上是一个群论概念，但是在 $\mathcal{OI}$ 中，其常见的形式为其在数论中的特例：

给定非负整数 $a, b, p$ ，要求最小的非负整数 $x$ ，满足 $a^x \equiv b \pmod{p}$ 。

接下来我们介绍两种用于求解数论中的离散对数的算法。

BSGS 算法

$\mathcal{BSGS}$ 全称 $\mathcal{Baby~Step, Giant~Step}$ ，又名 “大步小步” 算法，它实际上是基于一个类似分段打表的思想的。

实际上，根据底数 $a$ 和模数 $p$ 是否互质，这个算法又会有一些细节上的不同。

底数和模数互质

设 $q = \lfloor\sqrt{p}\rfloor, x = q \cdot r - s$ ，则方程可以化为：

$a^{q \cdot r - s} \equiv b \pmod{p}$

进一步有：

$(a^q)^r \equiv b\cdot a^s \pmod{p}$

因为 $s \in [0, q)$ ，所以我们可以考虑枚举 $s$ ，然后把所有 $b \cdot a^s \bmod p$ 用哈希表存储（如果值有重复，那么选择最大的 $s$ 即可，这样可以保证答案最小），再从小到大枚举 $r$ ，如果对于当前枚举到的 $r$ ， $(a^q)^r$ 可以在哈希表中查找到，那么此时的 $q \cdot r - s$ 就是最小的 $x$ 。

由于 $x$ 是 $\mathcal{O}(p)$ 级别的，因此 $r$ 同 $s$ 一样，也是 $\mathcal{O}(\sqrt{p})$ 级别的，因此算法的时间复杂度为 $\mathcal{O}(\sqrt{p})$ 。

底数和模数不互质

注意到此时上述的两个同余方程并非互为充要条件，因此我们需要先对方程做一些转化。

首先特判 $b = 1$ 的情况，此时原方程有解 $x = 0$ 。

设 $g = \gcd(a, p)$ ，那么根据取模运算的性质，原方程可以转化为：

$a^{x - 1}\cdot\frac{a}{g} \equiv \frac{b}{g} \pmod{\frac{p}{g}}$

显然，若某次转化后有 $g \nmid b$ ，那么原方程是无解的。

如果保证有解，那么我们可以考虑不断地转化这个方程，直到 $a \perp \dfrac{p}{\displaystyle \prod_{i = 1}^k g_i}$ 。

此时有：（设 $\displaystyle t = \prod_{i = 1}^k\dfrac{a}{g_i}$ ）

$\displaystyle a^{x - k} \equiv \dfrac{b}{\displaystyle t \cdot \prod_{i = 1}^k g_i} \equiv \frac{b}{a^k} ~(\bmod~\frac{p}{\displaystyle \prod_{i = 1}^k g_i})$

这个时候就可以用普通的 $\mathcal{BSGS}$ 求解了，最后的答案记得要加上 $k$ 。

当然，算法过程中还需要特判 $\displaystyle t = \frac{b}{\displaystyle \prod_{i = 1}^k g_i}$ 的情况，此时 $a$ 的次数为 $0$ ，答案即为 $k$ 。

Pohlig-Hellman 算法

这个算法一般用于求解 $p$ 为质数且 $p - 1$ 为 $\mathcal{Smooth~Number}$ （指质因数分解过后每个质因数均较小的数）的情况。但实际上，只要 $p$ 存在原根且 $\varphi(p)$ 为 $\mathcal{Smooth~Number}$ ，这个算法就是适用的。

下文默认 $p$ 为质数且且 $p - 1$ 为 $\mathcal{Smooth~Number}$ 。

首先考虑一个模 $p$ 的原根 $g$ ，设 $a \equiv g^{a'} \pmod{p},b \equiv g^{b'}\pmod{p}$ ，那么 $a^x \equiv b \pmod{p}$ 等价于 $g^{a'\cdot x} \equiv g^{b'}\pmod{p}$ ，也即 $a' \cdot x \equiv b' \pmod{p-1}$ 。

所以我们仅需求出 $a'$ 和 $b'$ ，便可以用 $\operatorname{exgcd}$ 求出 $x$ 的最小非负整数解。

现在我们聚焦于问题：对于方程 $g^x \equiv a\pmod{p}$ ，求出其最小非负整数解。

首先将 $p - 1$ 质因数分解，然后考虑对每一个 $u^v$ ，求出 $x \bmod{u^v}$ ，最后再用中国剩余定理合并。

我们将 $x$ 表示为 $u$ 进制，即 $\displaystyle x = (\sum_{i = 0}^{v - 1}r_i\cdot u^i) \bmod{u^v}$ ，那么问题在于如何求出所有的 $r_i$ 。

考虑构造方程 $(g^x)^{\frac{p - 1}{u^s}}\equiv a^{\frac{p - 1}{u^s}} \pmod{p^q}$ ，然后将里面的 $x$ 展开可得 $(g^{r_0 + r_1 \cdot u + r_2 \cdot u^2 +\cdots + r_{v - 1}\cdot u^{v - 1}})^{\frac{p - 1}{u^s}}\equiv a^{\frac{p - 1}{u^s}} \pmod{p^q}$ 。

若我们令 $s = 1$ ，那么原方程可以转化为 $g^{\frac{r_0 \cdot (p - 1)}{u}} \cdot g^{r_1\cdot ( p - 1)} \cdot g^{r_2\cdot u \cdot ( p - 1)} \cdot \cdots \cdot g^{r_{v - 1}\cdot u^{v - 2} \cdot (p - 1)}\equiv a^{\frac{p - 1}{u^s}} \pmod{p^q}$ ，根据费马小定理，我们还可以将其转化为 $g^{\frac{r_0 \cdot (p - 1)}{u}}\equiv a^{\frac{p - 1}{u^s}} \pmod{p^q}$ ，到了这一步，我们直接枚举 $r_0$ ，或是用 $\mathcal{BSGS}$ 算法求解 $r_0$ 均可。

接下来我们再令 $s = 2 \sim v$ ，分别求出 $r_{1\sim v - 1}$ 的值，即可求出 $x$ 的值。

显然时间复杂度和 $u$ 的大小有关，因此本算法适用于 $p - 1$ 为 $\mathcal{Smooth~Number}$ 的情况。

此处有一个小问题：为什么一定先把 $a^x \equiv b\pmod{p}$ 转化为 $a \equiv g^{a'} \pmod{p}$ 和 $b \equiv g^{b'}\pmod{p}$ 再求解呢？不能直接将其套用上文的算法求解吗？

实际上，在上文的算法中，我们利用到了 $\delta_p(g) = p - 1$ 这一条性质。如果 $\delta_p(g) < p - 1$ ，那么对于上文中的方程 $g^{\frac{r_0 \cdot (p - 1)}{u}}\equiv a^{\frac{p - 1}{u^s}} \pmod{p^q}$ ，若存在 $g^{\frac{p - 1}{u}}\equiv 1 \pmod{p^q}$ ，那么我们根本无法求出 $a_0$ 。

二次剩余

对于正整数 $n$ ，若 $x ^ 2 \equiv n \pmod{p}$ 有解，则称 $n$ 为模 $p$ 的二次剩余，否则为模 $p$ 的二次非剩余。

下文我们仅讨论 $p$ 为奇质数的情况。对于其它的情况，笔者并不了解相关算法，因此略过不谈。

注意当 $n = 0$ 时方程必然有解，此时我们认为 $0$ 是模 $p$ 的二次剩余，但下文提到 “二次剩余” 时，我们默认不考虑这种情况。

解的情况

我们先来讨论一下解的情况。

设方程有两个不等的解，分别为 $x_1,x_2$ ，则有 $x_1^2 \equiv x_2^2\pmod{p}$ ，移项后有 $(x_1 - x_2) \cdot (x_1 + x_2) \equiv 0\pmod{p}$ ，而由于前者必然不等于 $0$ ，所以 $x_1$ 和 $x_2$ 在模 $p$ 意义下必为相反数，并且由于 $p$ 为奇质数，因此它们是互不相等的。

也就是说，方程最多有两互异解，并且它们互为相反数。同时，由于我们不考虑 $n = 0$ 的情况，因此有解的方程一定有恰好两互为相反数的解。

进一步地，由于一对相反数仅会对应一个二次剩余，那么二次剩余的个数为 $\dfrac{p - 1}{2}$ ，其它的数即为二次非剩余，二者的数量相等。

Legendre 符号与 Euler 判别准则

$\mathcal{Legendre}$ 符号记作 $\displaystyle \left(\frac{a}{p}\right)$ 。若其等于 $1$ ，则说明 $a$ 为模 $p$ 的二次剩余；若其等于 $-1$ ，则说明 $a$ 为模 $p$ 的二次非剩余；若其等于 $0$ ，则说明 $a$ 在模 $p$ 意义下与 $0$ 同余。

值得一提的是， $\mathcal{Legendre}$ 符号还有更多的性质，比如它实际上是一个完全积性函数。感兴趣的读者可以自行查阅相关资料。

至于如何对 $\mathcal{Legendre}$ 符号求值，对于 $a \equiv 0 \pmod{p}$ 的情况我们无需在意，而对于剩下的情况，我们需要判别 $a$ 是否为模 $p$ 的二次剩余。

$\mathcal{Euler}$ 判别准则就是因此而生的：

当 $x ^ {\frac{p - 1}{2}} \equiv 1\pmod{p}$ 成立时， $x$ 为模 $p$ 的二次剩余，否则为二次非剩余，即 “ $x ^ {\frac{p - 1}{2}} \equiv 1\pmod{p}$ 成立” 是 “ $x$ 是模 $p$ 的二次剩余” 的充要条件。

证明：

先证必要性，设 $x \equiv y ^ 2 \pmod{p}$ ，由费马小定理可知 $x ^ {\frac{p - 1}{2}} \equiv y ^ {p - 1} \equiv 1 \pmod{p}$ 。

再证充分性，设 $g$ 为模 $p$ 的原根，那么存在 $i$ 使得 $x \equiv g^i \pmod{p}$ ，又因为 $x^{\frac{p - 1}{2}} \equiv g^{\frac{i\cdot (p - 1)}{2}} \equiv 1 \pmod{p}$ ，根据阶的性质， $\dfrac{i\cdot(p - 1)}{2}$ 一定是 $p - 1$ 的倍数，所以 $i$ 是偶数，令 $m = g^{\frac{i}{2}}$ 就有 $x \equiv m^2 \pmod{p}$ 。

综上，我们仅需判断 $x ^ {\frac{p - 1}{2}} \equiv 1\pmod{p}$ 是否成立，就可以知道 $x$ 是否为二次剩余。

Cipolla 算法

判定了 $n$ 是否为模 $p$ 的二次剩余后，我们怎样求出方程的解呢？

下面笔者将介绍 $\mathcal{Cipolla}$ 算法。

对于方程 $x^2 \equiv n \pmod{p}$ ，我们先随机出一个 $y$ 使得 $y ^ 2 - n$ 为二次非剩余，由于其个数与二次剩余相等，因此期望 $2$ 次就能找到。

然后我们类似虚数，将数域扩充，定义 $i ^ 2 \equiv y ^ 2 - n \pmod{p}$ 。

接下来给出结论： $x = (y + i) ^ {\frac{p + 1}{2}}$ 为原方程的一个解。

证明：

首先我们可以通过 $\mathcal{Lucas}$ 定理证明 $(y + i) ^ p \equiv y^p + i ^ p \pmod{p}$ 。

其次由费马小定理可知 $y ^ p \equiv y \pmod{p}$ 。

然后根据 $\mathcal{Euler}$ 判别准则和二次探测定理可知二次非剩余 $x$ 满足 $x^{\frac{p - 1}{2}} \equiv -1 \pmod{p}$ ，那么有 $(i ^ 2) ^ {\frac{p - 1}{2}} \equiv i ^ {p - 1} \equiv -1\pmod{p}$ ，即 $i ^ p \equiv -i \pmod{p}$ 。

所以：

$\begin{aligned} x ^ 2 &\equiv (y + i)^{p + 1} \\ & \equiv (y + i) ^ p \cdot (y + i) \\ & \equiv (y ^ p + i ^ p)\cdot (y + i) \\ & \equiv (y - i)\cdot (y + i) \\ & = y ^ 2 - i ^ 2 = n \end{aligned}$

我们求出 $x$ 之后，再通过上文中的结论求出另外一个解即可。

Miller-Rabin 质数判别算法

对于正整数 $n$ ，我们都知道存在时间复杂度为 $\mathcal{O}(\sqrt{n})$ 和 $\displaystyle\mathcal{O}(\sqrt{n}) - \mathcal{O}(\frac{\sqrt{n}}{\ln n})$ 的判别其是否为质数的算法，但是对于较大的 $n$ ，这个算法的时间复杂度还是比较难以接受的。

是否有复杂度更优秀的算法呢？

费马测试

我们不妨考虑质数的性质，其中有一条名为费马小定理：若 $p$ 是质数，则对于所有 $x \in [1, p)$ 都有 $x^{p-1} \equiv 1 \pmod{p}$ 。

那么，我们是否可以通过随机一些 $x$ ，判断上述条件是否对于这些 $x$ 均成立，从而判断 $p$ 是否为质数呢？

很遗憾，不行，原因是存在卡迈克尔数，对于这些数，我们无论如何都会将它们判别为质数。虽然它们的数量很少，但要预处理并特判它们也并不是简单的。

二次探测定理

我们不妨再来考虑质数的性质，其中有一条名为二次探测定理：若 $p$ 是质数，且有 $x^2 \equiv 1 \pmod{p}$ ，则 $x = 1$ 或 $x = p - 1$ 。

那么，我们是否可以通过将 $p - 1$ 分解为 $q \cdot 2^r$ 的形式，然后随机一些 $x$ ，依次判断是否有 $x^{q \cdot 2^{s + 1}} \equiv 1 \pmod{p}$ 且 $x^{q \cdot 2^s} \equiv 1~\texttt{or}~p -1 \pmod{p}$ ，最后再判断是否有 $x^{p -1}\equiv 1 \pmod{p}$ ，从而判断 $p$ 是否为质数呢？

很幸运，这是可以的，在 $x$ 随机的情况下，这个算法被证明有极高的正确率。事实上，参考这里和这里，我们可以通过调整测试数以得到确定性的判断结果，这样我们就无须担心算法的正确率，而仅需考虑一些需要特判的情况了。

最后，这个算法是可以做到单次测试 $\mathcal{O}(\log_2n)$ 的时间复杂度的，读者不妨自行思考如何实现。

需要代码实现的读者可以看这里。

Pollard-Rho 因数分解算法

对于正整数 $n$ ，我们还是都知道存在时间复杂度为 $\mathcal{O}(\sqrt{n})$ 、 $\mathcal{O}(n) - \mathcal{O}(\log_2n)$ 和 $\displaystyle \mathcal{O}(\sqrt{n}) - \mathcal{O}(\frac{\sqrt{n}}{\ln n})$ 的对其进行质因数分解的算法，但是对于较大的 $n$ ，这些算法的时间复杂度还是比较难以接受的。

是有有复杂度更优秀的算法呢？

确定方向

同 $\mathcal{Miller-Rabin}$ 算法，我们可以考虑在非确定性算法上寻找出路：我们可以通过随机化算法找出 $n$ 的一个非平凡因子 $m$ ，即满足 $m \ne 1~\texttt{or}~n$ 且 $m \mid n$ 。这样，若找出 $p$ 的复杂度为 $\mathcal{O}(f(n))$ ，那么算法的总时间复杂度为 $\mathcal{O}(f(n)\log_2n)$ 。

直接随机显然是不行的，我们考虑一种更为优秀的随机化算法。

生日悖论

我们来思考这样一个问题：不考虑出生年份，并假设每年均为 $365$ 天，那么当一个房间内有至少多少人时，有两个人生日相同的概率能够达到 $50\%$ ？

直觉告诉我们答案应该为 $183$ 人，但如果我们列出不等式，便可以得到一个十分反直觉的正确答案： $23$ 人。这也就是为什么该数学事实被称为一个悖论。

注意到 $23$ 是一个 $365$ 的平方根级别的数。实际上，根据生日悖论，我们可以导出一个推论：若我们在 $[1, n]$ 的范围内随机生成一个数并将其加入序列中，那么当序列中第一次出现两个相等的数时，序列的长度是期望 $\mathcal{O}(\sqrt{n})$ 级别的。

这个推论将在下文起到非常重要的作用。

构造伪随机函数

我们考虑构造这样一个伪随机函数： $f(x) = (x^2 + c)\bmod n$ ，其中 $c$ 是随机的常数。

这个函数有什么性质呢？

我们随机一个初值 $s$ ，然后观察序列 $s, f(s), f(f(s)), \cdots$ ，可以发现，这个函数会从某个位置开始进入一个循环。如果将有向图模型建出，可以发现其和 $\rho$ 的形状相类似，这便是这个算法的名字的由来。

设上述序列为 $x_i$ ，再设序列 $y_i$ ，且有 $y_i = x_i \bmod m$ （其中 $m$ 为 $n$ 的最小非平凡因子），可以发现 $y_{i +1} = (y_i^2 +c) \bmod m$ 仍然成立，也即 $y_i$ 的有向图模型同样和 $\rho$ 的形状相类似。

注意到若存在一对位置 $i, j$ ，满足 $x_i \not\equiv x_j \pmod{n}$ 和 $y_i \equiv y_j \pmod{m}$ ，则有 $n \nmid |x_i - x_j|$ 及 $m \mid |x_i - x_j|$ ，那么此时 $\gcd(n, |x_i - x_j|)$ 即为 $n$ 的一个非平凡因子。

而根据上文生日悖论的推论， $|i - j|$ （也即 $y_i$ 的 $\rho$ 形图的环长）期望是 $\mathcal{O}(\sqrt{m})\approx \mathcal{O}(n^{\frac{1}{4}})$ 级别的，那么，如果我们能够做到在 $\mathcal{O}(|i - j|)$ 的时间复杂度内找到 $n$ 的一个非平凡因子，我们就得到了一个复杂度十分优秀的算法。

算法流程

接下来我们考虑如何实现这个算法。

我们可以设置两个值 $u$ 和 $v$ ，分别从初值 $s$ 开始，每次令 $v$ 等于 $f(v)$ ，然后判断 $\gcd(n, |u - v|)$ 是否为一个非平凡因子。

然而这个算法存在两个问题：其一是 $u$ 可能不在 $y_i$ 的 $\rho$ 形图的环上，其二是有可能我们遍历了整个 $x_i$ 的 $\rho$ 形图的环都无法成功地找到一个非平凡因子。

这两个问题均可以用 $\mathcal{Floyd}$ 判环法来解决。具体地，我们还是设置两个值 $u$ 和 $v$ ，分别从初值 $s$ 开始，每次我们令 $u$ 等于 $f(u)$ ，令 $v$ 等于 $f(f(v))$ 。

这样，我们期望在 $\mathcal{O}(\sqrt{m})$ 次后遍历整个 $x_i$ 的 $\rho$ 形图的环，此时要么我们已经找到了一个非平凡因子，要么我们需要重新执行一遍上述算法。

这个算法的时间复杂度已经很优秀了，但其中包含了较多次数的 $\gcd$ 运算，这实际上会大大拖慢程序的运行速度。

根据相关性质，我们可以考虑把一部分的 $|u - v|$ 乘在一起，即令 $\displaystyle mul = (\prod |u - v|) \bmod n$ ，然后判断 $\gcd(n, mul)$ 是否为一个非平凡因子。具体地，我们可以考虑倍增 $|i - j|$ 的范围，然后将对应范围内的 $|u - v|$ 乘在一起，经实测当倍增到 $128$ 时停止倍增有最优的运行效率。

这样，我们可以把 $\gcd$ 运算的次数看做一个常数，那么算法的期望时间复杂度即为 $\mathcal{O}(\sqrt{m})$ 。

实际上， $\mathcal{Pollard-Rho}$ 算法的复杂度并没有得到十分严谨的证明，其原因是伪随机函数是否足够均匀还有待商榷。但不可否认的是，算法的运行效率确实非常高，这也让它能够很好地胜任大整数分解的任务。

需要代码实现的读者可以看这里。

积性函数的亚线性筛法

我们都知道埃氏筛的时间复杂度是近线性的，欧拉筛的时间复杂度是线性的，二者均可以筛出范围内的所有函数值。

但如果我们仅需求出函数值的前缀和，并且 $n$ 的范围能够达到 $10^9$ 甚至更高的级别，那么上述两种筛法就不太适用了。

接下来笔者将会介绍三种筛法，它们能够做到在亚线性的时间复杂度内求出积性函数的前缀和。

杜教筛

假设现在我们要求 $f(n)$ 的前缀和，即 $\displaystyle S(n) = \sum_{i = 1} ^n f(i)$ 。

我们考虑找出一个函数 $g(n)$ ，把 $f(n)$ 和 $g(n)$ 做狄利克雷卷积，设其前缀和为 $S'(n)$ ，那么有：

$\begin{aligned} S'(n) &= \sum_{i = 1}^n(f \cdot g)(i) \\ &= \sum_{i = 1} ^n \sum_{j | i} f(j) \cdot g(\frac{i}{j}) \\ &=\sum_{i = 1}^n g(i) \sum_{i = j}^{\lfloor\frac{n}{i}\rfloor} f(j) \\ &=\sum_{i = 1}^ng(i)\cdot S(\lfloor\frac{n}{i}\rfloor) \end{aligned}$

进一步地，有：

$g(1)\cdot S(n) = \sum_{i = 1}^ng(i)\cdot S(\lfloor\frac{n}{i}\rfloor) - \sum_{i = 2}^ng(i)\cdot S(\lfloor\frac{n}{i}\rfloor)$

即：

$S(n) = \frac{\displaystyle S'(n)- \sum_{i = 2}^ng(i)\cdot S(\lfloor\frac{n}{i}\rfloor)}{g(1)}$

所以，只要找出合适的 $g(n)$ ，就能够快速求出 $S(n)$ 了。

举两个例子：

对于 $\mu(n)$ 的前缀和，我们发现有 $\mu * 1 = \epsilon$ ，即设 $g(n) = 1$ ，代入公式有 $\displaystyle S(n) = 1 - \sum_{i = 2}^nS(\lfloor\frac{n}{i}\rfloor)$ 。

对于 $\varphi(n)$ 的前缀和，我们发现有 $\varphi \cdot 1 = \operatorname{id}$ ，即设 $g(n) = 1$ ，代入公式有 $\displaystyle S(n) = \frac{n(n + 1)}{2} - \sum_{i = 2}^n S(\lfloor\frac{n}{i}\rfloor)$ 。

这些递推式显然可以整除分块做。

实际上，我们可以先预处理出 $n^{\frac{2}{3}}$ 以内的 $S(n)$ ，对于剩下的 $S(n)$ 用记忆化搜索处理。可以证明，这样可以做到最优的时间复杂度 $\mathcal{O}(n^{\frac{2}{3}})$ 。

Min_25 筛

这种筛法对于 $f(n)$ 有一定的要求：当 $n$ 为质数时， $f(n)$ 是一个关于 $n$ 的度数较小的多项式。

首先考虑求质数的前缀和：

$\sum_{i = 1,i \in \mathcal{Prime}}^n f(i)$

我们进行这样一个转化：考虑当 $n$ 为质数时 $f(n)$ 的运算公式，然后把所有数都看作质数并套用这个运算公式。不妨设这个新的函数为 $h(n)$ 。

前面提到，当 $n$ 为质数时， $f(n)$ 是一个度数较小的多项式。所以，我们可以把 $h(n)$ 拆成若干项，分开来算，最后将它们合并在一起即可。为什么要将其拆成若干项？因为若我们不考虑符号和系数，那么多项式的每一项都是完全积性函数，这将为接下来的运算提供很大的帮助。

设 $g(n, j)$ 为所有质数和 $1\sim n$ 中最小质因子大于第 $j$ 个质数的合数的 $h(n)$ 的和，并设 $p_j$ 为第 $j$ 个质数。

$g(n, 0)$ 是很好求的，套用公式或用拉格朗日差值均可，那么问题在于如何转移。

转移的基本思路是按照 $j$ 从小往大转移。

若 $p_j > \sqrt{n}$ ，此时不会再有更多的合数对值造成影响，因此有 $g(n, j) = g(n, j - 1)$ 。

否则有 $p_j \le \sqrt{n}$ ，考虑需要被减去的贡献，显然这些数一定含有质因子 $p_j$ ，所以我们需要减去 $\displaystyle h(p_j)\cdot g(\lfloor\frac{n}{p_j}\rfloor,j−1)$ ，但 $\displaystyle g(\lfloor\frac{n}{p_j}\rfloor,j−1)$ 中包含了 $\displaystyle \sum_{i=1}^{j - 1}h(p_i)$ 这一部分小于 $p_j$ 的质数的函数值，它们实际上不应被减去（因为 $p_i\cdot p_j$ 的最小质因子为 $p_i$ ），所以再把它们加上就好了。

总结可得递推式：

$\begin{aligned} g(n, j) &= g(n, j - 1)~(p_j > \sqrt{n}) \\ g(n, j) &= g(n, j - 1) - h(p_j)\cdot (g(\lfloor\frac{n}{p_j}\rfloor,j - 1) - \sum_{i=1}^{j - 1}h(p_i))~(p_j \le \sqrt{n}) \end{aligned}$

求出来的 $g(n, |\mathcal{Prime}|)$ 即为 $1 \sim n$ 内所有质数的函数值之和。

接下来，类⽐ $g(n, j)$ ，设 $S(n, j)$ 表⽰最⼩质因⼦大于第 $j$ 个质数的所有数的 $f(n)$ 的和。

对于质数的部分，答案即为 $\displaystyle g(n, |\mathcal{Prime}|) - \sum_{i = 1}^{j} f(p_i)$ ，那么我们仅需考虑合数的部分。

枚举合数时，我们需要考虑所有最小质因子为 $q~(q > p_j)$ 的数 $r$ 对应的 $f(r)$ 的贡献，而显然合数 $r$ 的最小质因子是不超过 $\sqrt{r}$ 的，所以我们仅需考虑满足 $q^2 \le n$ 的部分即可。

对于每一个质因子 $q$ ，我们还需考虑其指数 $h$ 。由于 $S(n, j)$ 中实际上不包含 $f(1)$ 的贡献，那么我们需要先加入 $f(q^h)$ 这一部分贡献（注意这里不能重复加入质数的贡献）。对于剩下的部分，我们类似 $g(n, j)$ 的处理，利用积性函数的性质将 $f(q^h)$ 这一项提出来，再递归求解所有最小质因子大于 $q$ 的 $S(n, j)$ 即可。

总结可得递推式：

$S(n, j) = g(n, |\mathcal{Prime}|) - \sum_{i = 1}^{j} f(p_i) + \sum_{k = j + 1}^{p_k^2 \le n}\sum_{h = 1}^{p_k^{h} \le n} (f(p_k^{h + 1}) + f(p_k^{h}) \cdot S(\lfloor\frac{n}{p_k^h}\rfloor, k + 1))$

（注意当 $p_k^{h + 1}$ 大于 $n$ 时，我们认为 $f(p_k^{h + 1})$ 为 $0$ 。）

最后得到的 $S(n, 0)$ 即为答案，不过还需要加上 $f(1)$ 的贡献。

具体到代码实现上还有较多的细节，建议读者参考网上其它大佬的代码实现。

Powerful Number 筛

对正整数 $n$ 质因数分解为 $\displaystyle\sum_{p_i\in \mathcal{Prime}} p_i^{q_i}$ ，如果对于所有的 $i$ ，有 $q_i \ge 2$ ，则称 $n$ 为 $\mathcal{Powerful~Number}$ 。

实际上我们可以证明： $1 \sim n$ 内 $\mathcal{Powerful~Number}$ 的数量在 $\mathcal{O}(\sqrt{n})$ 级别。而这一条性质可以保证 $\mathcal{Powerful~Number}$ 筛的复杂度。

假设我们要求积性函数 $f(n)$ 的前缀和，那么我们需要找出一个函数 $g(n)$ ，满足其前缀和很好求，且两函数在质数处值相等，再找出一个积性函数 $h(n)$ 使得 $f(n) = h(n) \cdot g(n)$ 。

可以发现，当 $n$ 不为 $\mathcal{Powerful~Number}$ 时， $h(n)$ 取值为 $0$ 。

又有：

$\begin{aligned} \sum_{i = 1}^n f(i) &= \sum_{i = 1}^n (h \cdot g)(i) \\ &=\sum_{i = 1}^n\sum_{j | i} h(j)\cdot g(\frac{i}{j}) \\ &=\sum_{i = 1}^nh(i)\sum_{j = 1}^{\lfloor\frac{n}{i}\rfloor} g(j) \end{aligned}$

由于 $g(n)$ 的前缀和很好求，因此我们暴力地找出 $n$ 以内的 $\mathcal{Powerful~Number}$ ，算出对应的 $h(n)$ 之后再将它们相加即可。

本文没有讲解这些筛法的时间复杂度及其证明，对这些内容感兴趣的读者可以参考 $2018$ 年朱震霆的国家集训队论文《一些特殊的数论函数求和问题》。