Pollard-Rho算法证明

2020-12-12

昨天晚上简单地学习了一下Pollard-Rho算法，感觉很是神奇。同寝室的哲明宝宝真的太强了，虽然并没有接触过计算机竞赛，但是却对计算机竞赛数学掌握得炉火纯青，令我叹服。

借着那股热情，我和哲明宝宝讨论算法直到破晓。在哲明的帮助下，我觉得我已经能够比较详尽地给出Pollard-Rho算法的证明了，但是由于我现在表达能力严重退化，很多语言表述是很不严谨的，希望同学们谅解。

定理一：几何概型数学期望的证明

对于某次独立重复随机试验，实验成功的概率是 $p$，实验失败的概率是 $1-p$，倘若我们持续实验，直到第一次实验成功，试验次数的期望为 $\frac{1}{p}$。

这条定理几乎每个人在高中都学过，但是真正要是证明出来可能还是需要点大一的简单的数学才行，在这里给出一种可能比较麻烦的证法，这是我和室友哲明两个人的共同成果。

首先设 $E(p)$ 为单次实验成功概率为 $p$ 时，试验次数的数学期望，先给出一种很简单的证明方式。

$$E(p)=1\times p+(E(p)+1)\times (1-p)$$

解得 $E(p)=\frac{1}{p}$，轻松愉快。

这种方法很好理解，如果第一次实验就成功了，实验总次数为1。如果第一次实验没有成功，那么我们除了第一次实验外，就还需要进行 $E(p)$ 次实验。

此时我们给出更直观的做法，写出 $E(p)$ 的定义：

$$E(p)=p+2p(1-p)+3p(1-p{)}^2+\cdots$$

这是无穷项的和式，观察到指数与系数间存在差一的关系，想到对 $E(p)$ 进行积分。

$$\int (k+1)p(1-p{)}^k\text{d}p=?$$

可以采用换元，令 $q=1-p$，原始即为：

$$\int (k+1)(1-q)q^k\text{d}(1-q)=-q^{k+1}+\frac{k+1}{k+2}\cdot q^{k+1}+C$$

这个时候，就可以得到：

$$\int E(p)\text{d}p=\frac{1}{2}{p}^2+\sum _{k=1}^{\infty }\left(-(1-p{)}^{k+1}+\frac{k+1}{k+2}\cdot (1-p{)}^{k+2}\right)$$

将 $k=1$ 时的 $-(1-p{)}^{k+1}$ 从和式中提出。

$$E(p)=\frac{1}{2}{p}^2-(1-p{)}^2-\sum _{k=3}^{\infty }\frac{(1-p{)}^k}{k}$$

然后我们看到了惊人的一幕，我们看到了 $\frac{q^k}{k}$，这难道不正是$\ln (1+x)$在 0 处的麦克劳林展开式去掉前两项吗？

$$\ln (1+x)=x-\frac{1}{2}{x}^2+\frac{1}{3}{x}^3-\frac{1}{4}{x}^4+\cdots$$

因此，带入$x=-x$：

$$\ln (1-x)=-x-\frac{1}{2}{x}^2-\frac{1}{3}{x}^3-\frac{1}{4}{x}^4-\cdots$$

因此，带入$x=(1-p)$：

$$\ln (1-(1-p))=\ln (p)=-(1-p)-\frac{1}{2}(1-p{)}^2-\sum _{k=3}^{\infty }\frac{1}{k}(1-p{)}^k$$

解得：

$$\sum _{k=3}^{\infty }\frac{1}{k}(1-p{)}^k=\ln (p)+(1-p)+\frac{1}{2}(1-p{)}^2$$

因此：

$$\int E(p)\text{d}p=\frac{1}{2}{p}^2-(1-p{)}^2+\ln (p)+(1-p)+\frac{1}{2}(1-p{)}^2+C$$

对等式右侧求导可得：

$$E(p)=\frac{1}{p}$$

这个证明方式虽然看起来没有前面的那个证明简洁，但是前面的证明实际上并不严谨，为什么？因为 $E(p)$ 根据定义，本身就是一个无穷项的和式，只有当 $E(p)$ 极限存在的时候，我们才可以使用上述的证明方式，但是我们并没有证明极限存在，所以，第一种证明方式并不严谨。

假设一：随机数列假设

这一部分是我们全篇证明中最不美好的地方，因为这里我们给出了一个假设，而没有给出证明：

对于给定的正整数 $N,x_0,C$，构造数列：

$x_k=(x_{k-1}^2+C)\mathrm{m}\mathrm{o}\mathrm{d}N$

我们可以证明，数列 { x n } \{x_n\} {xn​} 一定存在周期，证明略。

还有一点，当 $x_0$ 和 $C$ 随机选择时，我们构造出的 { x n } \{x_n\} {xn​}在大多数情况下，在每个周期内，可以视为随机数列。这里的表述是十分模糊的，但是却很重要，后面的证明中会用到这个性质。

引理一：生日悖论

屋子里凑了一波乌合之众，然后它们把2月29号出生的人都踢了出去，好了，现在屋子里每种生日的人出现的概率都相同的。

假如说，现在，屋子里有366个人，那我我敢断言，物资中至少有2个人的生日在同一天。但是如果屋子里每个人的生日都是随机的，那么请问：如果我想让屋子中至少存在两个人生日相同这一事件发生的概率超过$\frac{1}{2}$，那么请问我们需要多少个人呢？如果您看过算法导论的话，您就会脱口而出：23个人就够了。

定义$P(N,M)$ ：

我有一个数列，数列的长度是 $N$，数列是一个随机数列，数列每一个位置上的数都是 $[0,M-1]$ 闭区间内的随机整数，那么该数列中至少存在两个相同的整数的概率我们就记为 $P(N,M)$。

显然：

$$P(N,M)=\frac{M^N-A_M^N}{M^N}$$

由于该数列符合古典概型，任意两种不同的数列方案，出现的概率相同，所以上式的正确性显然。

$M^N$ 表示数列的总方案数；$A_M^N$是一个排列，表示 $N$ 个数互不相同的方案数。

针对生日悖论，我们可以进行一些计算：

$$P(23,365)\approx 50.7\%$$

对于，一般情况，当 $M$ 很大的时候，我们总有：

$$P(\sqrt{M},M)\approx 1-\frac{1}{e}$$

这个东西的证明比较粗犷，用到了斯特林近似。

$$n!\approx \sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n$$

下面是我给出的非常不靠谱的推导。

$$P(N=\sqrt{M},M)=1-\frac{M!}{(M-N)!\cdot M^N}$$

记 $T=\frac{M!}{(M-N)!\cdot M^N}$。

$$T\approx \frac{\sqrt{M}\cdot {\left(\frac{M}{e}\right)}^M}{\sqrt{M-N}{\left(\frac{M-N}{e}\right)}^{M-N}\cdot M^N}$$

取对数，有：

$$T\approx e^{(M-N+\frac{1}{2})(\ln (M)-\ln (M-N))-N}$$

观察分子中的 $\ln (M)-\ln (M-N)$可以用拉格朗日中值定理进行代换。

$$\ln (M)-\ln (M-N)=N\cdot \frac{1}{\xi },\xi \in [M-N,M]$$

而当 $N=\sqrt{M}$ 时，就有：

$$\frac{1}{\sqrt{N}-1}=\frac{N}{M-N}\le \ln (M)-\ln (M-N)\le \frac{N}{M}=\frac{1}{\sqrt{M}}$$

因为 $M$ 很大，所以 $\frac{1}{\sqrt{M}}\approx \frac{1}{\sqrt{M}-1}$。

因此：

$$\ln (M)-\ln (M-N)\approx \frac{1}{\sqrt{M}}$$

因此：

$T\approx e^{(M-\sqrt{M}+\frac{1}{2})/\sqrt{M}-\sqrt{M}}=e^{\frac{1}{2\sqrt{M}}-1}\approx \frac{1}{e}$

所以：

$$P(\sqrt{M},M)\approx 1-\frac{1}{e}$$

引理二：随机数列的期望循环节长度为 $O(\sqrt{N})$

根据前文，伪随机数列的生成算法：

$$x_k=(x_{k-1}^2+C)\mathrm{m}\mathrm{o}\mathrm{d}N$$

不难发现：

如果出现了 $i\ne j$ 使得 $x_i=x_j$ 那么 一定有 $\forall k,x_{i+k}=x_{j+k}$，也就是说，循环节的长度一定是 $j-i$ 的因子。

因为 $P(\sqrt{N},N)\approx 1-\frac{1}{e}=\frac{e-1}{e}$。我们假设，我们使用了 $k$ 个长度为 $\sqrt{N}$ 的数列，使得出现 $i\ne j$ 使 $x_i=x_j$的情况，根据几何分布期望公式，我们能证明$E(k)=\frac{1}{P(\sqrt{N},N)}=\frac{e}{e-1}$。

也就是说，倘若我们有一个长度为 $L$ 的数列，使得 $L>$ 循环节的长度，那么 $L$ 的数学期望一定要 $\le \frac{e}{e-1}\sqrt{N}$。

这里的不等号是可以证明的，在此只给出比较形象的解释长度为 $L$ 随机数序列可以看作若干个长度为 $\sqrt{N}$的随机数序列拼接而成。$L>$ 循环节长度的充要条件是 $x_1\sim x_L$ 中找到至少两个元素相同，而 $\frac{e}{e-1}\sqrt{N}$是说有 $\frac{e}{e-1}$ 个长度为 $\sqrt{N}$的$x_n$的子串，这些子串中至少存在一个子串使得子串内包含两个相同的元素，这一条件显然比 $L$ 要满足的条件苛刻，因此 $E(L)\le \frac{e}{e-1}\sqrt{N}=O(\sqrt{N})$。

引理三：质因数分解的最坏情况

朴素质因数分解的时间复杂度为 $O(\sqrt{N})$ 这是因为，当 $N$恰好是一个质数的平方的时候，或者恰好是两个比较接近 $\sqrt{N}$ 的质数的乘积的时候，我们必须枚举 $2\sim \sqrt{N}$中的所有整数进行试除，而这种情况也正是质因数分解时间的瓶颈，时间复杂度瓶颈的证明需要在Pollard-Rho算法给出后给出。

定理二：Pollard-Rho 的算法思想

不妨假设 $N=p\times q,p\le q$ 且 $p,q$均为质数。

随机选取 $x_0$ 和 $C$，构造数列 $x_k=(x_{k-1}^2+C)\mathrm{m}\mathrm{o}\mathrm{d}N$。如果发现 { x k } \{x_k\} {xk​}是常数列，就重新选取一组 $x_0$ 和 $C$。

我们记 $b_k=x_k\mathrm{m}\mathrm{o}\mathrm{d}p$，根据引理二，我们能够得到数列 { b k } \{b_k\} {bk​} 的循环节的期望为 $O(\sqrt{q})$，而数列 { x k } \{x_k\} {xk​} 的循环节的长度为 $O(\sqrt{N})$。不出意外的话， { b k } \{b_k\} {bk​} 的循环节长度 $T_b$ 一定小于 { x k } \{x_k\} {xk​}的循环节长度 $T_x$，而且 $T_b$ 一定是 $T_x$的因子，也就是 $\exists k\in N_{+},k\times T_b=T_x$。

一种意外情况就是 $(x_0^2+C)\mathrm{m}\mathrm{o}\mathrm{d}N=x_0$，很尴尬，找到的数列 { x k } \{x_k\} {xk​}是常数列。这里是我证明的一个硬伤：我暂时不能证明：除了 $T_x=T_n=1$的情况外，是否还存在$T_x=T_b$的情况，我们暂且认为没有，谁要是证明了或者证伪了欢迎评论。

我们能够看到，$b_k$ 的循环可以看作一个“$\rho$”型的链表，当$k\ge l$ 时 $b_k$进入循环。

由于 $T_b<T_x,T_b|T_x$，那么我么一定存在一对整数 $1\le i,j\le T_x$ 使得 $T_b|(j-i)$ 也就是 $b_i=b_j$。

这个时候神奇的事情出现了：

$$x_i-x_j=(v_i\times p+b_i)-(v_j\times p+b_j)=(v_i-v_j)\times p+(b_i-b_j)=(v_i\times p+b_i)-(v_j\times p+b_j)=(v_i-v_j)\times p$$

其中数列 $v_k$ 的定义如下：

$$v_k=\lfloor \frac{x_k}{p}\rfloor$$

由于 $x_i\ne x_j,b_i=b_j$ 因此有 $v_i\ne v_j$，证明可以用反证法。

而 $Pollard-Rho$ 的质因数分解的原理就是：当 $k\ge l$ 时，如果我们能枚举$b_i$末端环上的每一个元素，那么肯定能找到一对 $i\ne j$ 使得$b_i=b_j$，即$x_i-x_j$非零且是 $p$ 的整数倍，此时 $\gcd (N,x_i-x_j)=p$，就这样我们就能找到 $p$。

不难发现，我们并不需要知道 $p$ 的值，并不需要知道 $b_i$ 的值，只需要通过递推数列 { x k } \{x_k\} {xk​}，枚举 $x_k$ 中的位置做差，就能计算出 $p$。

如何以 $O(\sqrt{p})=O(N^{\frac{1}{4}})$的时间复杂度求出一组 $b_i=b_j$呢？下述算法（python风格的伪代码）真的是让人拍案叫绝。

while true:
	x = randint(2, N-1)
	y = x # 其实比较点
	c = randint(2, N-1)
	i = 1
	s = 2
	while true:
		x = (x*x + c)%N
		if y == x: # 发现了T_x=T_b的情况
			break
		p = gcd(abs(x-y), N)
		if 1<p and p<N:
			print("find p!")
			return p
		i ++
		if i == s: # 重置起始比较点 y
			s *= 2
			y = x # 这个部分是算法令人拍案叫绝的地方

不难发现上述伪代码中的变量 $y$ 和变量 $x$ 都是我们前文中数列 $x_k$ 的某一项。首先，我们为什么要每迭代若干次重置一下 $y$ ？这是一个非常好的问题。

还是看图：

最开始的时候 $y=x_0$ 而 $b_0$实际上并不在环上，我们必须得保证 $y$ 和 $x$ 都在环上，才能保证能找到 $\gcd (abs(x-y),N)=p$，否则您可能循环一辈子也找不到 $p$。

其次，为什么我们每隔 $2^v,v\in N_{+}$次迭代，重置一下 $y$ ？这是为了均摊时间复杂度，使得我们迭代的总次数不超过 $O(l+T_b)=O(T_b)$（显然$l<T_b$，否则前 $l$ 个元素内部 $b_k$ 就已经进入循环了，与图示情况不符）。

而这种倍增的思路，会在当 $2^v$ 第一次大于等于 $T_b$的时候，找到 $p$。为什么？因为 $2^v\le p$，本次重置 $y$之后，下次重置 $y$值之前，足够 $x$值绕长度为$T_b$ 的 “$\rho$”环走一圈了。

此时，看看我们的 $x$ 走了多少步：

$$=1+2+4+8+\cdots +2^{\lfloor {\log }_2{T}_b\rfloor }+T_b=2^{\lfloor {\log }_2{T}_b\rfloor +1}-1+T_b=O(T_b)$$

而 $E(T_b)=\sqrt{q}=O(N^{\frac{1}{4}})$。

所以：$Pollard-Rho$ 的时间复杂度为 $O(N^{\frac{1}{4}})$，常数目测不小。

证明引理三

我们可以证明，倘若 $N$ 是质数或者 $N$是三个以上质数的乘积，时间复杂度一定不足 $O(N^{\frac{1}{4}})$。

首先，当 $N$ 是一个质数的时候，我们可以使用米勒拉宾素数测试以 $O(\log N)$ 的时间复杂度证明其为质数。

其次，当 $N$ 是三个质数的乘积的情况下，我们可以这样表述：$N=p_1{p}_2{p}_3$ 其中，$p_q,p_2,p_3$均为质数且 $p_1\le p_2\le p_3$。

我们用$Pollard-Rho$算法能在$O(\sqrt{p_1})$的时间复杂度找到质因子 $p_1$，此时由于 $p_1\le p_2\le p_3$，所以 $p_1\le N^{\frac{1}{3}}$，$\sqrt{p_1}=O(N^{\frac{1}{6}})$。

找到质因子 $p$ 后，我们还可以用 $O\left({\left(\frac{N}{p_1}\right)}^{\frac{1}{4}}\right)\le O\left(N^{\frac{1}{4}}\right)$ 的时间复杂度找到另外的两个质因子，不难证明 $O\left({\left(\frac{N}{p_1}\right)}^{\frac{1}{4}}\right)+O(N^{\frac{1}{6}})=O(N^{\frac{1}{4}})$。

综上所述，当 $N$ 的质因子个数更多的时候显然也能得到时间复杂度为$O(N^{\frac{1}{4}})$。