因式分解

Published: July 21, 2020 by Daltao

• Categories:
• algorithm 101

• Tags:

• Pollard-Rho大整数分解

Pollard-Rho大整数分解

对于任意正整数$n$，由算术基本定理可知，$n=\prod_{i=1}^kp_i^{c_i}$，其中$p_i$为互不相同的素数。

而因式分解是指，将一个合数分解为两个非平凡因子的乘积，比如$6=2\cdot 3$。

如何对合数$n$进行因式分解呢？

最简单的方式是枚举它可能的因子，由于$n$的最小因子一定在$2$到$\sqrt{n}$之间，因此可以在时间复杂度$O(\sqrt{n})$内完成。

Pollard's Rho算法在1975年由John Pollard发明，对于分解合数$n$，它的期望时间复杂度仅为$O(n^{\frac{1}{4}})$。

首先我们考虑建立一组序列$X=(x_1,x_2,\ldots)$，其中$x_{i+1}=x_i^2+c\pmod n$。由于是处于$n$的剩余类环中，因此序列中实际包含的不同值最多为$n$个。序列的图像应该如下：

由于图像酷似一个$\rho$，因此算法得名rho。

由生日悖论可知，在一年有$d$天的情况下，当班级里有$\sqrt{d}$个学生的时候，学生中期望至少有一对学生的生日相同，即发生碰撞。将生日悖论用在我们的序列中，可知我们序列中的不同值的数量约莫为$O(\sqrt{n})$。

设$m$为$n$的最小非平凡因子，易知$m\leq \sqrt{n}$。我们记$y_i=x_i\pmod m$，可以推导：

$$\begin{aligned} y_{i+1}&=x_{i+1}\pmod m\\ &=(x_i^2+c\pmod n)\pmod m\\ &=x_i^2+c\pmod m\\ &=(x_i\pmod m)^2+c\pmod m\\ &=y_i^2+c\pmod m \end{aligned}$$

因此我们得到了新的序列$Y=(y_1,y_2,\ldots)$，同样根据生日悖论，可知序列中不同值的数量约为$O(\sqrt{m})\leq O(\sqrt[4]{n})$。

由于第一次碰撞等概率发生在序列上，因此序列$X$中的环的入口期望为第$\sqrt{n}/2$个元素，而$Y$中的环的入口期望为第$\sqrt{m}/2$个元素。因此$X$中环的长度$C_x$期望为$\sqrt{n}/2$，同理$Y$中的环的长度$C_y$期望为$\sqrt{m}/2$。由于环的大小期望不同，那么就应该存在一对数$i,j$,满足$x_i\neq x_j \land y_i=y_j$。这意味着$n\nmid \left| x_i - x_j \right| \land m \mid \left| x_i - x_j \right|$。因此我们可以通过$gcd(n,\left|x_i - x_j \right|)$获得一个$n$的非平凡因子。

那么我们该如何获得$i$和$j$呢。有两种方法，一个是Floyd判圈算法，另外一种是倍增法。当Floyd判圈法记录两个下标$i$,$j$，每次循环$i=i+1$，$j=j+2$。当$C_y \mid (j-i)$时，我们就能找到对应的$i$和$j$。而倍增法则通过枚举可能的$Y$序列中环上的某个值和环的长度上限实现，首先猜想$y_1$落在环上且环的长度不大于$1$，如果第$i$次猜测失败，下一次就猜测$y_{2^i}$落在环上且环的长度不大于$2^i$。两种算法都可以在$O(\sqrt{m})$次循环内完成。

如果你挑选的$x_1$和$c$不好的话，可能会导致最终没法找到$n$的任意因子。那么你可以在合适的时机，用新的随机数替代$x_1$和$c$重新进行启发式查找。

注意，这个算法是建立在$n$是合数的前提下的，如果$n$是素数，那么这个算法就永远无法找到因子而陷入死循环中。所以在执行Pollard's Rho算法之前，需要先测试$n$是否是合数，这里可以采用随机算法Miller Rabin进行，这里不过多讲。

并且这里$O(\sqrt{m})$是指的是循环次数，而不是时间复杂度，由于每次循环我们都会调用gcd，因此真实的时间复杂度为$O(\sqrt{m}\log_2 n)$。

int pollard_rho(int x, int c, int n)
{
    int xi = x;
    int xj = x;
    int j = 2;
    int i = 1;
    while(i < n)
    {
        i++;
        xi = (xi * xi + c) % n;
        int g = gcd(n, abs(xi - xj));
        if(g != 1 && g != n)
        {
            return g;
        }
        if(i == j)
        {
            j = j << 1;
            xj = xi;
        }
    }
    return -1;
}

上面的是倍增法的代码。