理论基础

  阅读这篇文章，需要有双指针查环的算法基础。可以看我写的两篇博文Floyd链表查环算法和Brent链表查环算法。为了辅助理解，还可以先读读我写过的线性同余伪随机数算法。
  Pollard $\rho$算法是使用同余算法产生一个$[0,n)$之间的伪随机数，这里的n是我们要进行分解的数字，不过不是这个算法线性同余，可以说是二次同余，因为用的函数是$f(x)=(x^2+c)modn$。
  根据生日悖论birthday paradox，这个随机数序列组成的链表一定会成环。我举个例子，$c=1,x_1=3,p=72$。会生成如下序列：

  这个图形就长得像字母$\rho$，所以叫$\rho$算法。这个环的检测使用Brent或Floyd链表查环算法就行了。但是对于更大的数进行质因数分解，这个环会非常大，比如这个数2206637，如果对2206637取模，形成的环一个屏幕是装不下的。所以对这个序列中的值还可以继续对质因子取模。因为$6961=317\times 6961$，我们以$c=1,x_1=2,p=317,n=2206637$为例子，画一下对原序列再进行一次取模运算的效，这里我选择对小的质因子317进行取模：

  在上图中冒号前的是原值，冒号后的是对p再次取模后的值。169那个节点有个括号，表示取模前两个值是不同的，但是取模后都是169。那存在环意味着什么呢？这和质因数分解有什么关系呢？假设相遇的这个点和相遇前的一个点，在序列里的顺序（索引）分别为$s,t$，以1开始。比如上图，$s=6,t=12$。再设序列为 { x    m o d    p } \{x \;mod\; p\} {xmodp}，注意这里的x是取模前的值，$x_s$和$x_t$是对$p$取模同余的，所以有：
$$\begin{gathered} x_t\equiv x_{s}modp \\ |x_t-x_s|=kp,k\in \mathbb{Z} \end{gathered}$$
  用上面的例子代进去更容易理解，上例中$x_6=1166412,x_{12}=1264682$（用前一项计算的结果），有：
$$|x_6-x_{12}|=98270=310\ast 317=310p$$
  因为p是我们需要求出的一个质因数，所以有：
$$n=pq$$
  如果形成了环，那么意味着，$|x_t-x_s|=kp$，而$kp$和$n=pq$的最大公约数就是p的倍数。到这一步，意味着我们检验成环的条件变了，变成了求和n的最大公约数了。如果最大公约数不是1，那么这个数字就是n的一个因子。那我们就找到了一个，然后n再除于这个最大公约数，这样递归下去就完成了质因数分解。用公式总结就是：
$$\begin{gathered} x_t\equiv x_{s}modp \\ |x_t-x_s|=kp,k\in \mathbb{Z} \\ gcd(|x_t-x_s|,n)=p \end{gathered}$$
  然后我们再使用查环算法进行查环，代码就出来了。

通用代码

def gcd(m, n):
    if m < n:
        m, n = n, m
    while m % n != 0:
        m, n = n, m % n
    return n


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


def factorization(n, fac_fun):
    factors = []
    stack = [n]
    while len(stack) > 0:
        x = stack.pop()
        if x == 2:
            factors.insert(0, x)
            continue
        p, q = fac_fun(x) if x & 1 == 1 else (2, x // 2)
        if p == 1:
            factors.insert(0, q)
        else:
            stack.append(p)
            stack.append(q)
    return factors
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

Floyd查环代码实现

def floyd(n):
    c = 1
    x = 2
    slow = x
    fast = x
    while True:
        slow = f(slow, c, n)
        fast = f(f(fast, c, n), c, n)

        diff = abs(fast - slow)
        if diff != 0:
            g = gcd(diff, n)
            if g > 1:
                return g, n // g
        else:
            return 1, n

if __name__ == '__main__':
    n = 2206637
    print(factorization(n, rho))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

  测试结果如下：

[317, 6961]
1

  注意代码中的细节，如果是对质数进行Pollard $\rho$算法分解，相差会得到0。

Brent查环代码实现

def brent(n):
    c = 1
    x = 2
    slow = x
    fast = f(x, c, n)
    power = 1
    steps = 1
    while True:

        diff = abs(fast - slow)
        if diff != 0:
            g = gcd(diff, n)
            if g > 1:
                return g, n // g
        else:
            return 1, n

        if steps == power:
            power <<= 1
            steps = 0
            fast = slow

        slow = f(slow, x, n)
        steps += 1


if __name__ == '__main__':
    n = 2206637
    print(factorization(n, brent))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

  测试结果如下：

[317, 6961]
1