Pollard’s rho算法

引入

给定一正整数$N\in {\mathbb{N}}^{\ast }$，试快速找到它的一个因数。
很久很久以前，我们曾学过试除法来解决这个问题。很容易想到因数是成对称分布的：即$N$的所有因数可以被分成两块：$[1,\sqrt{N}]$和$[\sqrt{N},N]$。这个很容易想清楚，我们只需要把区间$[1,\sqrt{N}]$扫一遍就可以了，此时试除法的时间复杂度为$O(\sqrt{N})$。
对于$N\ge 10^{18}$的数据，这个算法运行起来无疑是非常糟糕的。我们希望有更快的算法。对于这样的算法，一个比较好的想法是去设计一个随机程序，随便猜一个因数。如果你运气好，这个程序的时间复杂度下限为$O(1)$。但对于一个$N\ge 10^{18}$的数据，这个算法给出答案的概率是$\frac{1}{1000000000000000000}$。当然，如果在$[1,\sqrt{N}]$里面猜，成功的可能性会更大。那么，有没有更好的改进算法来提高我们猜的准确率呢？那就是Pollard’s rho算法。

算法流程

Pollard’s rho算法的重要思想就是：最大公约数一定是某个数的约数。也就是说，$\forall k\in {\mathbb{N}}^{\ast }$，$gcd(k,n)|n$。只要选适当的k使得$gcd(k,n)>1$就可以求得一个约数$gcd(k,n)$。这样的k数量还是蛮多的：k有若干个质因子，而每个质因子的倍数都是可行的。

如果这样的k我们随机取的话，就体现不出Pollard’s rho算法的精妙之处了。我们不妨考虑构造一个伪随机数序列，然后取相邻的两项的差来求gcd。为了生成一串优秀的随机数，Pollard设计了这样一个函数：$f(x)=(x^2+c)\mathrm{mod}N$，其中$c$是一个随机的常数。

我们随便取一个$x_1$，令$x_2=f(x_1)$，$x_3=f(x_2)$，……，$x_i=f(x_{i-1})$。在一定的范围内，这个数列是基本随机的，可以取相邻两项作差求gcd。

这样构造的数列会有一些性质：
（1）相邻两数的差会包含前面所有相邻两数差的乘积，这些乘积中可能会包含n的因子。考虑下面这个式子：
$$x_i-x_{i-1}=f(x_{i-1})-f(x_{i-2})={x_{i-1}}^2-{x_{i-2}}^2=(x_{i-1}+x_{i-2})(x_{i-1}-x_{i-2})(\mathrm{mod}n)$$
在等式最右边$x_{i-1}$和$x_{i-2}$又可以继续分解，直到$x_1$。所以会有上面那个性质。

（2）序列中距离为k的两项的差，一定为前面任意距离为k的两项的差的倍数。这个性质算是上面那个性质的推广，当$k=1$时，其实和性质（1）是一样的。在序列中取$x_i,x_j$，假设$i>j$，则有：
$$x_i-x_j=f(x_{i-1})-f(x_{j-1})={x_{i-1}}^2-{x_{j-1}}^2=(x_{i-1}+x_{j-1})(x_{i-1}-x_{j-1})(\mathrm{mod}N)$$
其中等式最右边$x_{i-1}$和$x_{j-1}$又可以继续分解。

（3）序列中任意两数的差，也一定可以转换为相邻两个数的差的倍数。在序列中取$x_i,x_j$，假设$i>j$，则有：
$$x_i-x_j=x_i-x_{i-1}+x_{i-1}-\cdots -x_j+x_j-x_{j-1}=b\ast (x_j-x_{j-1})(\mathrm{mod}N)$$
根据性质（1），相邻两数的差会包含前面所有相邻两数差的乘积。所以，$x_i-x_{i-1}$一直到$x_{j+1}-x_j$会包含$x_j-x_{j-1}$，将所有的无关项提取出来形成一个b就得到了上面等式的最右边。

根据前面的分析，序列有一些特殊的性质，即下标距离为k的两个数之差，是前面的下标距离为k的某两个数之差的倍数。所以，对于距离为k的两个数的差，我们只需要检测最后那一对即可。这样，每一个距离都只需要检测一次。

为了判断环的存在，可以用一个简单的Floyd判圈算法，也就是"龟兔赛跑"。假设乌龟为$t$，兔子为$r$，初始时$t=r=1$。假设兔子的速度是乌龟的二倍。过了时间$i$后，$t=i,r=2i$。此时两者得到的数列值$x_t=x_i,x_r=x_{2i}$。假设环的长度为$c$，在环内恒有：$x_i=x_{i+c}$。 如果龟兔"相遇"，此时有：$x_r=x_t$，也就是$x_i=x_{2i}=x_{i+kc}$。此时两者路径之差正好是环长度的整数倍。

出现这种乌龟被套圈的情况，再继续下去只会一直重复，所以需要通过调整函数$f(x)$里常数项的值来继续“碰运气”。

这样一来，我们得到了一套基于Floyd判圈算法的Pollard Rho 算法。Python实现如下

import random
from math import *
from Crypto.Util.number import isPrime


def f(x, c, n):
    return (x * x + c) % n


def Pollard_rho(n):
    if isPrime(n):
        return n
    while True:
        c = random.randint(1, n - 1)
        t = f(0, c, n)
        r = f(t, c, n)
        while t != r:
            d = gcd(abs(t - r), n)
            if d > 1:
                return d
            t = f(t, c, n)
            r = f(r, c, n)
            r = f(r, c, n)


n = int(input())
fac = Pollard_rho(n)
print(fac)

椭圆曲线上的Pollard’s rho算法

设椭圆曲线为$G$，$P,Q$为椭圆曲线上两点，欲求出满足$dP=Q$的$d$。

选择哈希函数将$G$分成尺寸大致相同的三部分$S_1,S_2,S_3$，其中$O\notin S_2$；定义一个随机步数的迭代函数$f$：
$$R_{i+1}=f(R_i)=\{\begin{array}{cc}Q+R_i,& R_i\in S_1\\ 2R_i,& R_i\in S_2\\ P+R_i,& R_i\in S_3\end{array}$$
令$R_i=a_i\ast P+b_i\ast Q$，则：
$$a_{i+1}=\{\begin{array}{cc}{a}_i,& R_i\in S_1\\ 2a_i\mathrm{mod}N,& R_i\in S_2\\ a_i+1,& R_i\in S_3\end{array}$$
$$b_{i+1}=\{\begin{array}{cc}{b}_i+1,& R_i\in S_1\\ 2b_i\mathrm{mod}N,& R_i\in S_2\\ b_i,& R_i\in S_3\end{array}$$
初始化参数$R_0=P,a_0=1,b_0=0$，产生配对$(R_i,R_{2i})$ ，直到对某个$m$有$R_m=R_{2m}$，此时有：
$$R_m=a_{m}P+b_{m}Q$$
$$R_{2m}=a_{2m}P+b_{2m}Q$$
所以有$d=\frac{a_{2m}-a_m}{b_m-b_{2m}}\mathrm{mod}N$