N^{1/5} polylog N 因数分解 - Harvey 算法 简化版

· · 算法·理论

算法挺简单，实现起来略复杂。

前情提要

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

Harvey 算法 Slides Paper 两处的算法描述是不同的，不妨看本文吧。

Harvey 算法简明易懂的理论分析 https://www.luogu.com.cn/article/ocaz51pd 需要边思考边看或者先看过前文或者本文才好懂。

算法描述

简单的数学部分

令 N = pq、N^{2/5} < p \le N^{1/2}、x = up + vq、y=up - vq、知 x^2 - y^2 = 4uvN。选定 (u, v) 求解 x 即可求解 p。

如 \frac{u}{v} 是 v \le V = \frac{p}{(4N)^{2/5}} 的从某一侧最接近 \frac{q}{p} 的分数，有 |y| < (4N)^{2/5} = m。这是根据 Dirichlet 逼近定理：存在 v\le V 满足 |\frac{q}{p}-\frac{u}{v}|<\frac{1}{vV}。对于这个 (u, v) 求解 x 时，知道 x^2 \in [4uvn, 4uvn + m^2)。

现在对于未知的 (u, v)，令 z = uv \le up/m = (vq + y)/m = (N v/p + y) / m < N/m/m + 1，即 z \le Z = \frac{N}{m^2} = (4N)^{1/5} / 4。

对于每个可能的取值 z、知要求解的整数 x 可能在 [\sqrt{4zN}, \sqrt{4zN + m^2}) 之中，该区间长度小于等于 m^2\cdot\frac{1}{2\sqrt{4zN}}。当 z = \frac{N}{m^2} 时，区间长度小于等于 m^2\cdot\frac{1}{2\sqrt{4zN}} = \frac{m^3}{4N} = (4N)^{1/5}。

所有区间总长度小于等于 \sum_{z=1}^{Z}{m^2}\cdot\frac{1}{2\sqrt{4zN}} \approx m^2\sqrt\frac{Z}{4N} = \frac{m}{2}。考虑到 z = uv 的不同的 (u, v) 要对区间长度重复计数，区间总长度小于等于 \sum_{1 \le v \le u,\ uv \le Z}{m^2}\cdot\frac{1}{2\sqrt{4uvN}} = \sum_{v=1}^{\sqrt{Z}}\frac{m^2}{\sqrt{4vN}}\sum_{u=v}^{Z/v}\frac{1}{2\sqrt{u}} \approx \frac{m}{2}\sum_{v=1}^{\sqrt{Z}}\frac{1}{\sqrt{vZ}}\cdot(\sqrt{Z/v}-\sqrt{v}) = \frac{m}{2}\sum_{v=1}^{\sqrt{Z}}(\frac{1}{v}-\frac{1}{\sqrt{Z}})\approx \frac{m}{2} \ln Z。

算法

特判、初始化

首先假定 p > m = (4N)^{2/5}。否则使用 O(N^{1/4} \mathrm{\ polylog\ } N) 实际上是 O(\sqrt{p} \mathrm{\ polylog\ } N) = O(N^{1/5} \mathrm{\ polylog\ } N) 的算法因数分解。也可以先跑本算法，如果求不出解，再跑 O(\sqrt{p} \mathrm{\ polylog\ } N) 的算法。

选取 \alpha（方法见理论分析最后一段）使得 \mathrm{ord}_N(\alpha) > d（d 的值见下方第二段）、对于 k \le d 枚举 \gcd(\alpha^k \bmod N - 1, N) 判断是否为 p。不成立则知道 \mathrm{ord}_p(\alpha) > d。

详情

枚举整数 z \le Z = (4N)^{1/5} / 4。同时对每个区间 [\lceil\sqrt{4zN}\rceil, \sqrt{4zN + m^2})，按照下面的方法求解 x 或 p。

选取 d = [(4N)^{1/5}]、m \approx d^2、令 x_0 = \lceil\sqrt{4zN}\rceil、x_i = x_0 + id。形象地说，是对每个 [\lceil\sqrt{4zN}\rceil, \sqrt{4zN + m^2}) 区间砍成长度为 d 的段，最后一段长度不满 d 无所谓。

求解 $x$ 相当于求解 $\alpha^{j} = \alpha^{u+vN - id - x_0} = \alpha^{t(u,v,i)} \mod p$。左边为 $d$ 个不同数值、右边为 $\frac{m}{2d}\ln Z + Z \ln Z = 0.75 d \ln Z = F$ 个不同的 $t(u,v,i)$ 数值。详细解释一下 $F$ 的计算。第一项代表所有的区间，按照长度 $d$ 分段，大约是几段。$\ln Z$ 因子代表同一个 $z$ 对应了平均 $\ln Z$ 组 $(u, v)$。加法第二项代表每个区间的最后一段可能是不满长度 $d$ 的，因此对于每一段再加 1，一共是增加区间总量。 找出所有右边的 $\alpha^{t(u,v,i)} \bmod N$ 值，在 $\alpha^k \bmod N$ 中去重，丢掉重复数值后，称为 $t'(u,v,i)$，令多项式 $f(x) = \Pi_{t'(u,v,i)}(x-t'(u,v,i))\bmod N$。对于去重丢掉的 $\alpha^{t(u,v,i)} = \alpha^j$，使用 $p = \frac{x\pm\sqrt{x^2-4zn}}{2u}$ 计算出 $p$ 判断是否整数，其中 $x = x_0 + id + j$。 运用离散傅立叶变换类似的 [Chirp Z 变换]( https://oi-wiki.org/math/poly/czt/) 同时计算 $f(\alpha), \ldots, f(\alpha^d)$。之后计算 $\gcd(f(\alpha^j), N)$，如果不是 1，可能是 $p$ 或者 0。如果是 0，对于 $f(\alpha^j)$ 的每个乘数项计算 $\gcd(\alpha^j-t'(u,v,i), N)$，必有一个不是 1，而是 $p$。 因为 $f$ 的表达式是乘法形式，实际上要先计算出多项式系数，问题描述为知道它的根，求系数，进行分治和 DFT 解决。$f$ 的阶 $F$ 大约和 $d$ 只相差一个很小的 $\log N$ 因子。 时间复杂度是 $O(d \mathrm{\ polylog\ } N)$ 即 $O(N^{1/5} \mathrm{\ polylog\ } N)$。$\log N$ 因子次数是 3 吧。