模意义下高次方程

Prologue

数论科技越来越多我该怎么办？数论科技越来越多我该怎么办？数论科技越来越多我该怎么办？

本文主要将介绍模意义下的高次方程求解，包括 N 次剩余、离散对数、幂塔等。

都是基础知识，没有数论科技

Part 1. 开根号！！1

OI 中常见的是二次剩余和比较 general 的 N 次剩余。但其实也有三次剩余的快速求解方法。

1.1 二次剩余

求解 x^2\equiv a\pmod p。先只考虑 p 为奇素数的情况。

二次剩余的存在性

Theorem（Euler 判别准则）：对于 p\nmid a，有：

进行一个意识流的证明：

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}：

a^{(p-1)/2}=(g^v)^{(p-1)/2}=(g^{(p-1)/2})^v

显然有 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} 进行二项式展开：

\begin{aligned} (r+Y)^{p+1}=&\ \sum_{i=0}^{p+1}\binom{p+1}ir^iY^{p+1-i}\\ =&\ r^{p+1}+Y^{p+1}+(p+1)(r^pY+rY^p)\\ =&\ r^2-Y^2+(p+1)(rY-rY)\\ =&\ r^2-(r^2-a) =a \end{aligned}

接着要证明这样算出来的数中的虚部，即 \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}。

1. 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} 成立。直接转移即可。

2. 1. $b\equiv 0\pmod {p^\alpha}

   容易发现 a 必然为 p^{\lceil\frac\alpha n\rceil} 的倍数（包括 0），枚举即可。

   2. b\perp p

   找出 p^\alpha 的原根 g，令 g^\lambda\equiv x，g^\mu\equiv b，那么有 \lambda n\equiv \mu\pmod{\varphi(p^\alpha)}。exgcd 求解线性方程即可。

   3. 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 算法

令 $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](https://www.luogu.com.cn/problem/P4195) ### 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 离散对数](https://loj.ac/p/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}$。即： $$ p^{\alpha-\sigma-1}\log_gh=p^{\alpha-\sigma-1}\sum_{i=0}^{\sigma}x_ip^i $$ 移项后把这个式子扔到指数上可得： $$ \left(h\times g^{-\sum_{i=0}^{\sigma-1}x_ip^i}\right)^{p^{\alpha-\sigma-1}}=\left(g^{p^{\alpha-1}}\right)^{x_\sigma} $$ 由于 $(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](http://acm.hdu.edu.cn/showproblem.php?pid=6632) ## 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(\log m)$ 的。这个数列更严谨的上下界可见[OEIS](https://oeis.org/A003434)。 Exercise：[Luogu P4139 上帝与集合的正确用法](https://www.luogu.com.cn/problem/P4139) [Luogu P3747 [六省联考 2017] 相逢是问候](https://www.luogu.com.cn/problem/P3747) #### 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(\log p)$ 的递归上界。于是我们只需特判 $a=1,2$ 和 $a=b$ 的情况，否则直接按无限层的幂塔计算即可。 $n>3$ 时类似 $n=3$ 的做法，特判一些 corner cases 即可。 Exercise：[Luogu P6736 「Wdsr-2」白泽教育](https://www.luogu.com.cn/problem/P6736) ## References Part 1: [1](https://blog.csdn.net/zxyoi_dreamer/article/details/85195819) [2](https://www.sciencedirect.com/science/article/pii/S0893965902000319) [3](https://www.luogu.com.cn/blog/Leasier/solution-P5668) Part 2: [1](https://www.luogu.com.cn/blog/Leasier/index-calculus) [2](https://www.sciencedirect.com/science/article/pii/S0747717199902791) EOF