N^{1/5} polylog N 因数分解

· · 算法·理论

尽管容易理解，没有 OI 价值吧…… 简化版算法直接关掉本文看这里：https://www.luogu.com.cn/article/040zdn9a

前情提要

https://www.luogu.com.cn/article/9jfrrk1k Lehman 算法（1974）O(N^{1/3}) 实现

Pollard O(N^{1/4}) 算法基础上用 Lehman 算法的 O(N^{1/3}) 的算式（注意顺序）魔改的https://web.maths.unsw.edu.au/~davidharvey/talks/magma2023.pdf O(N^{1/5}) Harvey 算法。（O(N^{1/4}) 算法指的不是随机的这个Pollard-Pho）

看上面的文章时，发现很巧妙地避开了最佳逼近分数来实现，只是求解费马分解方程，因此对这一系列因数分解算法的很感兴趣。06 年做集训队作业的时候，研究 Farey 序列给出了一个 O(\log n) 求解任意实数附近最佳逼近分数的算法，这里用得上。注：（一类分数问题的研究.pdf 可能哪里找得到吧，e.g. https://blog.csdn.net/Code92007/article/details/135053808 但这里的不完整，没包括实数部分，实数到分数的方法只是先找一个平方规模的分母分子作为这个实数的代替求解，之后可以证明就是解）

Lehman 算法基本细节

[前情提要](https://www.luogu.com.cn/article/9jfrrk1k )里的初等观点，没有常数因子，避开了带来核心进展并且非平凡的最佳分数逼近，因此也失去了继续优化的空间。 这里只补充说明 [Slides](https://web.maths.unsw.edu.au/~davidharvey/talks/magma2023.pdf) 之中不显然的地方。因为懒得看论文原文，就自己推导了下。 ### 主方程、算法解释 稍稍给它改头换面，比较好从[前情提要](https://www.luogu.com.cn/article/9jfrrk1k )的理解跳过来。 $vN/x+ux = vq_0+up_0 + O(1)\ \forall x, s.t.\ |x-p_0| \le p_0/N^{1/3}, q_0 = N/p_0

意思是，考虑若干个 p_0，把 p \in [N^{1/3},N^{1/2}] 这个大区间，用 \forall x, s.t.\ |x-p_0| \le p_0/N^{1/3} 覆盖了。每个区间后端点对前端点的比值为 \frac{1+N^{1/3}}{N^{1/3}-1}，因此总区间数是\frac{\log\frac{N^{1/2}}{N^{1/3}}}{\log \frac{1+N^{1/3}}{N^{1/3}-1}} = O(N^{1/3}\log N)。

对于每个区间，令 $y = vq_0+up_0 + O(1)$ 进行枚举（只有 -1 到 +2 范围内的三个整数），求解 $p (=x) = \frac{y \pm \sqrt{y^2 - 4uvN}}{2u}$，隐含了根式首先是整数。跟[前情提要](https://www.luogu.com.cn/article/9jfrrk1k )中的求解方程 $x^2-y^2=4uvN$ 有同样的实质，只是方法差异很大。相当于对于每组 $(u, v)$，通过对应的 $p_0$ 直接确定了整数 $(x, x^2-4uvN)$，而不是小范围枚举 $x$。 意思是，$x^2-4uvN$ 不是在一个小范围内自由的，而是由 $p_0$ 或者说某种程度上由 $(u, v)$ 所确定的。这是因为 $x^2-y^2=4uvN$ 这个方程丢失了主方程中除法和整性所包含的值域信息。这个方程丢失信息的主要原因是只用了 $\frac{u}{v}$ 逼近在 $uv$ 小范围的存在性，而没考虑进每个 $(u,v)$ 锁定了一个 $p$ 的可行范围。 ### 主方程的详细说明 $vN/x+ux = vq_0(p_0/x)+up_0 (x/p_0) = (vq_0+up_0)+vq_0/x(p_0-x)+u(x-p_0)=(vq_0+up_0)+(x-p_0)/p_0(up_0-vq_0)+(x-p_0)^2vq_0/p_0x = (vq_0+up_0) + \frac{x-p_0}{p_0}\cdot\frac{p_0}{V}\cdot vV(\frac{u}{v}-\frac{N}{p_0^2})+\frac{(x-p_0)^2}{xp_0}N\frac{v}{V}\frac{V}{p_0}$。其中 $\frac{p_0}{V}$ 是选定的参数，$|\frac{p_0}{x-p_0}|$ 的下界是另一个参数。观察第二项和第三项，知两个参数都取 $N^{1/3}$ 时，均为 $O(1)$。 另一种参数选取，分别为 $N^{2/5}$、$N^{1/5}$ 时，均为 $O(N^{1/5})$。 这里注意一下用的是狄利克雷逼近定理来证明第二项的乘数因子 $O(1)$ 的。首先它是鸽笼原理证明的存在性定理。而正好它的 Farey 序列下一项 $a/b$ 的距离是 $1/(vb) \ge 1/(vV)$ 因此下一项跳过了被逼近的实数，因此只可能是单侧最优的，故“某一个”单侧最优的逼近的差值也满足狄利克雷逼近定理的界。 至于最好能做到 $1/5$ 次方，是因为第二项是 $a-b$ 次、第三项是 $1-2b-a$ 次，而时间复杂度中，有 $b$ 次个 $p_0$、而界之中是 $a-b$ 次，决定了 DFT 同时计算的多项式值的个数。故知 $\max(a-b, b)$ 最小为 $1/5 = [2(a-b)+3b]/5 = [1+(a-b)-(1-2b-a)]/5$。 ## $O(N^{1/5})$ Harvey 算法思想 因为假设了某个 $p_0$ 确定的区间存在 $p$，知道 $ vN/p+up - (vq_0+up_0) = O(N^{1/5})$， 对 $\alpha^{vN/p+up - \lceil vq_0+up_0\rceil} = \alpha^k \mod p$ 应用费马小定理，即为 $\alpha^{vN+u - \lceil vq_0+up_0\rceil} = \alpha^k \mod p$，左边次方上面没有未知数了。 左侧有最多 $O(N^{1/5}\log N)$ 个数值，令 $f(x) = \Pi_{i(p_i, u, v, q_i\mathrm{\ s.t. }\forall k\ \alpha^{vN+u - \lceil vq_i+up_i\rceil} \ne \alpha^k \bmod N)}(x-\alpha^{vN+u-\lceil vq_i+up_i\rceil}) \bmod N)$ 是 $O(N^{1/5}\log N)$ 次多项式。知道 $f(\alpha^k)$ 存在一个 $p$ 的倍数项。要计算 $f(\alpha^{-N^{1/5}}), \ldots, f(\alpha^{2N^{1/5}})$，之后对于每一个 $f(\alpha^k)$ 与 $N$ 求 gcd。 计算多个多项式的数值是 [Chirp Z 变换]( https://oi-wiki.org/math/poly/czt/)。 这部分是 Slides 所说的，但我看 [paper](https://arxiv.org/pdf/2010.05450) 采取的是完全不同的做法。不理它为好，这里只说明 Slides 的做法。 ### 因数分解的严谨性 在求出 gcd 之后，如果是 $p$，得到解了，如果是 $N$ （即 0），则逐个检查 $f(\alpha^k)$ 的因子和 $N$ 的 gcd，如果是 $p$，得到解了。如果是 $N$ 就会有麻烦。 前面知道存在 $(p_i, u, v)$，使得 $vN/p+up - \lceil vq_i+up_i\rceil = k$，而 $k = O(N^{1/5})$，又有假设，$\mathrm{ord}_N(\alpha) > O(N^{1/5})$，这里就不写精确的常数因子了，这里满足 $\mathrm{ord}_N(\alpha)$ 对应的 $> O(N^{1/5})$ 对应了所有可能 $k$ 中最大值。 因为可能 ，$vN+u - \lceil vq_j+up_j\rceil = k + l\cdot\mathrm{ord}_N(\alpha)$，这里 $p$ 并不在 $p_j$ 对应的区间之内。好的性质是，如果有这个情况，只可能对应一个 $k$。 对于这种情况，首先 $\{\alpha^k \bmod N\}$ 和 $\{\alpha^{vN+u-\lceil vq_i+up_i\rceil} \bmod N\}$ 都是预先计算的。先直接进行判重。知道对于 $p_j$，最多只有一个对应的 $k$，那么只需要对 $vN/p+up - \lceil vq_j+up_j\rceil = k$ 使用 $p = \frac{y \pm \sqrt{(k + \lceil vq_j+up_j\rceil)^2 - 4uvN}}{2u}$ 求解判定是否为整数。不论计算结果如何，这个 $(p_j, k)$ 这对关系就已经被排除掉了。我们最多排除掉 $O(N^{1/5} \log N)$ 对这样的关系，不影响总的时间复杂度。 现在考虑，如果对于 $(p_i, u, v)$，使得 $vN/p+up - \lceil vq_i+up_i\rceil = k$，并且 $vN+u - \lceil vq_i+up_i\rceil \ne k + j\cdot\mathrm{ord}_N(\alpha)$ 也就是我们希望 $p_i$ 应当可以解出 $p$ 的情况。考虑到是否存在 $k'$ 使得 $vN+u - \lceil vq_i+up_i\rceil = k' + l\cdot\mathrm{ord}_N(\alpha) = k \mod \mathrm{ord}_p(\alpha)$，如果存在，知 $\mathrm{ord}_p(\alpha) | k' - k = O(N^{1/5})$，即 $vN/p+up - \lceil vq_i+up_i\rceil = k = k' + c\cdot\mathrm{ord}_p(\alpha)$，不可以因为 $k'$ 的存在排除掉 $p_i$，同时知道 $1 < \gcd(\alpha^{k} - \alpha^{vN+u-\lceil vq_i+up_i\rceil}, N) < N$。 先从 $f(x)$ 中除去所有的存在 $\alpha^{k'} = \alpha^{vN+u-\lceil vq_i+up_i\rceil} \mod N$ 的所有 $p_i$ 有关的因子进行求解。如果找不到解，说明 $p$ 在某个这样的 $p_i$ 的区间之内。 对于这些 $\{p_i\}$，我们希望寻找 $k$ 使得 $1 < \gcd(\alpha^{k-k'} - \alpha^{vN+u-\lceil vq_i+up_i\rceil-k'}, N) < N$，即 $1 < \gcd(\alpha^{k-k'} - 1, N) < N$，这是平凡的，并且可以作为算法启动前的特判，即在相应的 $O(N^{1/5})$ 范围内朴素地枚举 $i$ 检验 $\gcd(\alpha^i - 1, N)$ 判段。 $\mathrm{ord}_p(\alpha) = O(N^{1/5})$ 时，可以直接得到因子 $p$。 因此算法在对 $f(x)$ 的乘法因子中排除某些区间后成立。 ### 其他细节补充 由于 $\frac{p_0}{V} = N^{2/5}$ 是算法选定的参数，隐含 $p > N^{2/5}$。如果刚刚求不出解，即落入此情况，使用 $O(N^{1/4})$ 算法先进行因数分解。注意到该算法的时间复杂度实际上是 $O(\sqrt{p}\log N) = O({N^{1/5}\log N})$。 ### $\alpha$ 的选取方法 要求 $\mathrm{ord}_N(\alpha) \ge D = k^2$。 对于某个 $x$，判断 $\{x^i \bmod N\} \cap \{x^{-jk} \bmod N\}$ 是否为空，其中 $1 < i, j < k$。空集说明 $\mathrm{ord}_N(x) \ge D$。否则知道出它的阶是 $\mathrm{ord}_N(x) = i+jk$。再选择 $a$，如果 $a^{\mathrm{ord}_N(x)} \ne 1$，求出阶，那么 $a$ 的阶是 $x$ 的阶的倍数，或者 $ab$ 的阶是两个阶的最小公倍数，否则，$a$ 没有用处。总之，新的阶至少是之前的两倍。 按照这个方法，令 $\alpha = p_1^{b_1} p_2^{b_2}\ldots p_{\log D}^{b_{\log D}} \ (b_i \in {0, 1})$，就是一个解。我们知道，根据算法的过程有 $\mathrm{ord}_N(p_i) ｜ \mathrm{ord}_N(\alpha) = M$，就知道了对于任意的 $2^{\log D}$ 种 $(a_i \in \{0, 1\})$ 的组合 $(p_1^{a_1} p_2^{a_2}\ldots p_{\log D}^{a_{\log D}})^M = 1 \mod q$。当 $p_1p_2\ldots p_{\log D} < q$ 时，这些数字互不相等，故此时多项式 $x^M - 1 = 0 \mod q$、$q \ge \sqrt{N}$ 有至少 $D$ 个不同的 1 之外的根，说明 $M > D$。不需要特地找素数，朴素枚举、没有新的素因子跳过就行了。 算法可以随时在 $\mathrm{ord}_N(\alpha) \ge D$ 时停止。我们这里的 $D$ 相对于 $q$ 足够小，上一段之中的 $p_1p_2\ldots p_{\log D} < q$ 基本上是成立的，因为左边小于 $4^{p_{\log D}} = O(N^{2/5} \log^2 N)$，右边至少 $\sqrt{N}$。其实还有更强的证明方式。 我们已经知道，在本算法范围内，上面的方法是可以找到 $\mathrm{ord}_N(\alpha) \ge D$ 的。在 $D$ 更大的情况，最坏情况是这样的：在总时间复杂度为 $O(\sqrt{D} \log^2 D)$ 的寻找后，如果找不到 $\alpha$，至少有 $\sqrt{D} \le M < D$，判断是否存在某个素因子 $p$ 使得 $\mathrm{ord}_p(\alpha)$ 更小，如果是可以通过 $\gcd(\alpha^i - 1, N)$ 得到 $p$。如果不是，知道 $N$ 的所有因子都是 $aM+1$，因此 $1 < \gcd(N, (M+1)(2M+1)\ldots([\sqrt{N}/M]M+1)) < N$，需要用到另一种质因数分解方法，大致可以参考以下文章稍作改动： "使用快速阶乘算法求出 B 阶下降幂"做因数分解 https://www.luogu.com/article/jyha8mh2、 https://www.luogu.com/article/0o2xs4uf （快速阶乘算法） 、 [快速点值平移](https://oi-wiki.org/math/poly/shift/#lagrange-%E6%8F%92%E5%80%BC%E5%85%AC%E5%BC%8F%E6%B3%95)。