模意义下高次方程：从开门到入门
Prologue
数论科技越来越多我该怎么办？数论科技越来越多我该怎么办？数论科技越来越多我该怎么办？
本文主要将介绍模意义下的高次方程求解，包括 \(N\) 次剩余、离散对数、幂塔等。
都是基础知识，没有数论科技
Part 1. 开根号！！1
OI 中常见的是二次剩余和比较 general 的 \(N\) 次剩余。但其实也有三次剩余的快速求解方法。
1.1 二次剩余
求解 \(x^2\equiv a\pmod p\)。先只考虑 \(p\) 为奇素数的情况。
二次剩余的存在性
Theorem（Euler 判别准则）：对于 \(p\nmid a\)，有：
\[a^{(p-1)/2}\bmod p=\left(\frac ap\right)=\begin{cases}1& a\text{ is a quadratic residue modulo }p\\-1 & \text{otherwise}\end{cases} \]其中 \(\left(\frac ap\right)\) 为 Legendre 符号。
没有卵用
进行一个意识流的证明：
Lemma：\(a\) 为模 \(p\) 意义下的二次剩余的充要条件为：令 \(g\) 为 \(p\) 的原根，则 \(\log_ga\) 为偶数。
略证：令 \(x=g^u\)，\(a=g^v\)，则有 \(2u\equiv v\pmod {(p-1)}\)。由于 \(2u\) 和 \(p-1\) 均为偶数，可得 \(v\) 为偶数。反之亦然。
考虑拆一下 \(a^{(p-1)/2}\)：
显然有 \(g^{p-1}=1\)，则 \(g^{(p-1)/2}=\pm1\)。由原根的定义可知 \(g^{(p-1)/2}=-1\)。然后讨论一下 \(v\) 的奇偶性就出来了。
由 Lagrange 定理，\(a^{(p-1)/2}\equiv1\pmod p\) 有 \(\frac{p-1}2\) 个解。因此模 \(p\) 意义下的二次剩余和非二次剩余个数均为 \(\frac{p-1}2\)。
Cipolla Algorithm
求解模奇素数意义下的二次剩余。
考虑定义一种类虚数。这种方法是求解 \(N\) 次剩余中的一种重要思考方法，下文的三次剩余求解也用到了这种思想。
考虑找到一个 \(r\)，使得 \(r^2-a\) 为模 \(p\) 意义下的非二次剩余。令 \(Y^2\equiv r^2-a\pmod p\)，并定义我们的类虚数为 \(\alpha+\beta Y\)。因为一个数为模 \(p\) 意义下二次剩余的概率为 \(\frac12\)，因此我们只需随机约 \(2\) 次便可找到一个 \(r\)。
算法流程很简单：定义一个类虚数 \(r+Y\)（即 \(\alpha=r\)，\(\beta=1\)），然后快速幂求解 \(\sqrt a\equiv(r+Y)^{(p+1)/2}\pmod p\)。
证明：首先要证 \((r+Y)^{p+1}=a\)。
Lemma：\(Y^{p-1}=-1\)。
略证：\(Y^{p-1}=(r^2-a)^{(p-1)/2}=-1\)。
考虑到模 \(p\) 意义下二项式系数中的阶乘可以直接被除掉，我们对 \((r+Y)^{p+1}\) 进行二项式展开：
接着要证明这样算出来的数中的虚部，即 \(\beta=0\)。考虑设 \((u+vY)^2\equiv a\pmod p,v\not=0\)。则有 \(u^2+v^2(r^2-a)+a=-2uvY\)。左边没有 \(Y\)，因此有 \(uv=0\)，\(u=0\)，即 \((vY)^2=a\)。如此可得 \(Y=v^{-1}\sqrt a\)，即 \(Y^2\) 为二次剩余，与条件矛盾。于是我们的算法是正确的。
时间复杂度就是计算快速幂的复杂度，为 \(O(\log p)\)。
不要尝试去碰任意模数二次剩余，会变得不幸。
Exercise：Luogu P5491 【模板】二次剩余
1.2 三次剩余
求解 \(x^3\equiv a\pmod p\)。仅考虑 \(p\) 为奇素数的情况。任意模数我不会啊，直接做 N 次剩余吧
首先特判一下 \(a\in\{0,1\}\text{ or }p=3\) 的情况，直接让 \(x=a\) 即可。
接着特判一下 \(p=3k+2\)，由费马小定理，直接令 \(x=a^{(2p-1)/3}\)。
然后是 \(p=3k+1\)。注意到 \(x^3\equiv1\pmod p\) 有非平凡根 \(\epsilon=\frac{-1+\sqrt{-3}}2\) 以及 \(\epsilon^2\)，因此有 \(\frac13\) 的数有三个模立方根，有 \(\frac23\) 的数没有立方根。
考虑令 \([a/p]=a^{(p-1)/3}\bmod p\)（下文证明会用到这个符号）。显然该式子的取值为 \(\{1,\epsilon,\epsilon^2\}\)。与 Euler 判别准则类似，当且仅当 \(a^{(p-1)/3}=1\) 时 \(a\) 为模 \(p\) 意义下的三次剩余。读者不妨模仿 Euler 判别准则的证明来证明这个结论。
Peralta Method
假设 \(a\) 为模 \(p\) 意义下的三次剩余。与 Cipolla Algorithm 类似，定义环为：令 \(Y^3=a\)，一个环可以表示为 \(z=\alpha+\beta Y+\gamma Y^2\)。由于 \(Y\) 对应 \(\mathbb Z_p\) 中的一个数，\(z\) 显然满足费马小定理，即 \(z^{p-1}\equiv 1\pmod p\)。
假设有一个 \(z\)，使得 \(z^{(p-1)/3}=\beta Y\)（\(\alpha=\gamma=0\)），那么显然有 \(x\) 的一个解为 \(\beta^{-1}\)（等式两边立方一下即可）。
那么问题变为了如何选取这样的 \(z\)。类似于 Cipolla Algorithm（对，又是类似），考虑随机选取 \(z\)。下面会进行一个意识流的证明：每次选取出满足条件的 \(z\) 的概率为 \(\frac 19\)。
考虑一个函数 \(\varphi\)：令 \(x_0\) 为 \(x^3\equiv a\) 的一个解，\(\varphi(z)=(\alpha+\beta x_0+\gamma x_0^2,\alpha+\beta x_0\epsilon+\gamma x_0^2\epsilon^2,\alpha+\beta x_0\epsilon^2+\gamma x_0^2\epsilon)\)，即 \(Y\) 取遍 \(x^3\equiv a\) 的三个解。那么 \(\varphi\) 为环的集合到 \(\mathbb Z_p\times\mathbb Z_p\times \mathbb Z_p\)（即一个模 \(p\) 意义下的整数三元组）的一个双射。手玩一下可以发现 \(\varphi(a+b)=\varphi(a)+\varphi(b)\)，\(\varphi(ab)=\varphi(a)\varphi(b)\)（三元组间的运算为对应位置上的数进行相应的运算）。
令 \(\varphi(z)=(r,s,t)\)，若 \(z^{(p-1)/3}=\beta Y\)，那么 \((r^{(p-1)/3},s^{(p-1)/3},t^{(p-1)/3})=(\beta x_0,\beta x_0\epsilon,\beta x_0\epsilon^2)\)。那么对于每个 \(i\in\{1,\epsilon,\epsilon^2\}\)，有 \(\beta=x_0^{-1}i\)，\([r/p]=i\)，\([s/p]=i\epsilon\)，\([t/p]=i\epsilon^2\)。对于 \(r\)（\(s,t\) 同理），由于规定了 \([r/p]\) 的取值，\(r\) 一共有 \(\frac{p-1}3\) 种取值。那么 \(\operatorname{Pr}(z^{(p-1)/3}=\beta Y)=\frac{3((p-1)/3)^3}{(p-1)^3}=\frac19\)。
论文上说这种算法的复杂度为 \(O(\log^3 p)\)，但是我怎么看都是 \(O(\log p)\) 的，希望有大佬能教一下。
读者不妨考虑一下如何将该算法扩展到求解二次剩余或五次剩余等。（实际上，该算法求解五次剩余及以上时常数和复杂度均较大，没啥用）
Exercise：Loj #175 模立方根
1.3 \(N\) 次剩余
求解 \(x^n\equiv b\pmod m\)。
考虑把 \(m\) 质因数分解，对其所有素数幂的因数求解，然后 CRT 合并。于是问题转化为求解 \(x^n\equiv b\pmod {p^\alpha}\)。
- \(p=2\)
考虑增量法构造。假设求出了 \(x^n\equiv b\pmod {2^{\alpha-1}}\) 在模 \(2^{\alpha-1}\) 意义下的解集 \(S\)，那么对于其中每个元素 \(x_0\)，在模 \(2^\alpha\) 意义下必定仅有两个解 \(x_0,x_0+2^{\alpha-1}\) 成立。直接转移即可。
- \(p\) 为奇素数
-
\(b\equiv 0\pmod {p^\alpha}\)
容易发现 \(a\) 必然为 \(p^{\lceil\frac\alpha n\rceil}\) 的倍数（包括 \(0\)），枚举即可。
-
\(b\perp p\)
找出 \(p^\alpha\) 的原根 \(g\)，令 \(g^\lambda\equiv x\)，\(g^\mu\equiv b\)，那么有 \(\lambda n\equiv \mu\pmod{\varphi(p^\alpha)}\)。exgcd 求解线性方程即可。
-
\(b\equiv 0\pmod{p}\)
若 \(x=\lambda p^e\)，\(b=\mu p^k\)，则有 \(\lambda^n\equiv \mu\pmod{p^{\alpha-k}}\)，且 \(en=k\)。于是转化成了情况 2。
需要注意多解的情况。
Exercise：Luogu P5668 【模板】N 次剩余 Luogu P8457 「SWTR-8」幂塔方程
Part 2. 取对数！！1
OI 中较常用 BSGS 算法和 exBSGS 算法。当然还有科技
2.1 BSGS & exBSGS
根号分治&预处理后根号分治。
BSGS 算法
求解 \(a^x\equiv b\pmod n\)，\(n\perp a\)。
令 \(m=\lceil\sqrt n\rceil\)，\(x=um-v\)。于是有 \((a^m)^u\equiv ba^v\pmod n\)。考虑预处理出 \(ba^v\) 对应的 \(v\) 值，然后枚举所有可能的 \((a^m)^u\) 值，寻找是否有 \(v\) 与之对应即可。复杂度 \(O(\sqrt p\log p)\)。
Exercise：Luogu P3846 [TJOI2007] 可爱的质数/【模板】BSGS Luogu P2485 [SDOI2011]计算器
exBSGS 算法
\(n\) 与 \(a\) 不互质时，发现 \(a^{um-v}\equiv b\Leftrightarrow (a^m)^u\equiv ba^v\) 的转化可能不成立。考虑把 \(n\) 与 \(a\) 的公因数给除掉。
令 \(d=\gcd(a,n)\)。$d\nmid b $ 时显然无解，否则令 \(n\mapsto n/d\)，\(b\mapsto b/d\)。假设如此迭代 \(k\) 次后 \(a\perp n\)，那么答案即为 \(\operatorname{BSGS}(a,b,n)+k\)，其中 \(\operatorname{BSGS}(a,b,n)\) 为使用 BSGS 算法算出的 \(a^x\equiv b\pmod n\) 的解。
Exercise：Luogu P4195 【模板】扩展 BSGS/exBSGS
2.2 Index Calculus 算法
求 \(a^x\equiv b\pmod p\)，\(p\) 为素数。
一句话概括：利用指数函数将原根转化为 smooth number 来计算出一些质数的前缀离散对数，再将要求离散对数的数转化为 smooth number 进行计算。
Smooth Number：称一个数为 \(N-\text{smooth}\)，当且仅当该数的所有质因子均 \(\le n\)。
考虑设定一个阈值 \(N\)，筛出 \(N\) 以内的所有质数。考虑求出 \(\log_g p\) 的值。考虑随机一个 \(x\)，若 \(g^x\) 为 \(N-\text{smooth}\)，那么将其质因数分解为 \(\prod p_i^{\alpha_i}\) 后取 \(\log\) 可得 \(\sum\alpha_i\log p_i\equiv x\)，在模 \(\varphi(p)\) 意义下高斯消元即可。注意到每行元素只有 \(O(\log)\) 个，可以用一下稀疏矩阵高斯消元的科技。
预处理后考虑解类似于 \(g^x\equiv b\pmod p\) 的方程。考虑随机 \(y\)，若 \(g^yb\) 是 \(N-\text{smooth}\) 的，那么将其质因数分解后取 \(\log\) 可得 \(x\equiv \sum\alpha_i\log p_i-y\)。
然后考虑 general 的情况 \(a^x\equiv b\pmod p\)。先解出 \(g^y\equiv a\) 以及 \(g^z\equiv b\)，那么有 \(xy\equiv z\pmod{\varphi(p)}\)。exgcd 求解同余方程即可。
论文指出当 \(N=\exp(\frac{\sqrt{\ln p\ln\ln p}}2+1)\) 时算法复杂度较优。复杂度证明要用到 quadratic sieve 相关理论，不会。
Exercise：Loj #6542 离散对数
2.3 Pohlig-Hellman 算法
求 \(a^x\equiv b\pmod p\)，\(p\) 为素数且 \(p-1\) 的质因数分解较简单。
类似于 Index Calculus 算法，将原问题转化为 \(g^x\equiv h\pmod p\)，即 \(\log_gh\equiv x\pmod{(p-1)}\)。
考虑将 \(p-1\) 质因数分解，于是问题转化为求解 \(\log_gh\equiv x\pmod {p^\alpha}\) 后 CRT 合并。把 \(x\) 写成 \(p\) 进制的形式，即 \(x=x_0+x_1p+x_2p^2+\dots+x_{\alpha-1}p^{\alpha-1}\)。考虑从前往后对 \(x_i\) 依次求解。
假设当前我们已经求解出了前 \(\sigma\) 项，即 \(x_0,x_1,\dots,x_{\sigma-1}\)（初始时 \(\sigma=0\)），要求 \(x_\sigma\)。发现如果在等式左右两边同时乘上一个 \(p^{\alpha-\sigma-1}\)，就能消去 \(x_{\sigma+1},\dots,x_{\alpha-1}\)。即：
移项后把这个式子扔到指数上可得：
由于 \((x_0x_1x_2\dots x_{\alpha-1})_p\) 是 \(x\) 的 \(p\) 进制表示，我们有 \(0\le x_\sigma<p\)，因此可以暴力枚举 \(x_\sigma\) 并快速幂检验。
算法的复杂度取决于 \(p-1\) 的最大质因子 \(\textit{maxp}\) 与质因子数量 \(\omega(p-1)\)。一般情况下的复杂度大概是 \(O(\textit{maxp}\log^2p\ \omega(p-1))\)。当 \(p\) 取某些特殊值（如 \(998244353=2^{23}\times7\times17+1\)，\(65537=2^{16}+1\)）时，算法复杂度约为 \(O(\log^2p)\)。
Exercise：HDU 6632 discrete logarithm problem
Part 3. 递归！！1
3.1 幂塔
求 \(a^{a^{a^{\dots}}}\bmod m\)。
考虑令原式等于 \(f(a,m)\)。由扩展欧拉定理可得 \(f(a,m)=a^{f(a,\varphi(m))+\varphi(m)}\bmod m\)。递归计算即可。
考虑这种算法递归执行次数的上界是多少，即至多进行多少次迭代能使得 \(\varphi(\varphi(\dots\varphi(m)))=1\)。对于一个奇数 \(m\)，\(\varphi(m)\) 显然为偶数（将其质因数分解后，必有一个奇素因子 \(p^k\) 对答案有 \(\varphi(p^k)\) 的贡献）；对于偶数 \(m\)，有 \(\varphi(m)\le \frac m2\)，因为每次取 \(\varphi\) 必定至少使 \(m\) 的一个素因子 \(2\) 变为 \(1\)。如此可知迭代次数是 \(O(\lg m)\) 的。这个数列更严谨的上下界可见OEIS。
Exercise：Luogu P4139 上帝与集合的正确用法 Luogu P3747 [六省联考 2017] 相逢是问候
3.2 高德纳箭号
对于 \(a,b\in \mathbb N\)，\(n\in\mathbb N^*\)，高德纳箭号定义为：
\[a\uparrow^nb=\begin{cases}1&b=0\\a^b&n=1\\a\uparrow^{n-1}(a\uparrow^n(b-1)) &\text{otherwise}\end{cases} \]
\(n=1\) 时即为指数运算，不再赘述。
\(n=2\) 时手玩一下发现 \(a\uparrow^2b=a\uparrow\uparrow b=a^{a^{a^{\dots}}}\)（共 \(b\) 个 \(a\)）。这时只需套用上面的幂塔算法，并且在递归时判断层数即可。
\(n=3\) 时，发现 \(a\uparrow^3b=a^{a^{a^{\dots}}}\)（共 \(a\uparrow^3(b-1)\) 个 \(a\)）。这个 \(a\) 的数量极大，极易超过 \(O(\lg p)\) 的递归上界。于是我们只需特判 \(a=1,2\) 和 \(a=b\) 的情况，否则直接按无限层的幂塔计算即可。
\(n>3\) 时类似 \(n=3\) 的做法，特判一些 corner cases 即可。
Exercise：Luogu P6736 「Wdsr-2」白泽教育
References
EOF
