Pollard-Rho 学习笔记

​
Pollard-Rho 是一种期望 $O(n^{\frac{1}{4}})$ 时间分解质因数的算法。

Pollard-Rho 用到了 Miller-Rabin，Miller-Rabin 是一种 $O(\log n)$ 时间判素数的方法。

首先来讲 Miller-Rabin，它用到了一种叫二次探测的东西和费马小定理。

Miller-Rabin

二次探测

有这样一个方程 $x^2\equiv 1(\mathrm{mod}p)$，$x$ 的根分别为 $1$ 和 $p-1$。

证明：

原方程即为 $(x-1)(x+1)\equiv 0(\mathrm{mod}p)$。

于是 $x=1$ 或 $x=-1$。

也就是 $x=1$ 或 $x=p-1$。

我们可以利用这个性质来判断 $x$ 是否为质数。

费马小定理

$x^{p-1}\equiv 1(\mathrm{mod}p)$

这个不用多说。

快速乘

因为要判断素数的数和要分解质因数的数都很大，乘一次就超出了 long long 范围，所以要用快速乘。

模板：

ll qpow(ll a, ll b, ll p) {
    return (a * b - (ll)((long double)a / p * b) * p + p) % p;
}

原理是余数等于被除数减商乘除数。

Miller-Rabin 算法

我们选取了一个质数 $p$，首先对 $x$ 判断费马小定理是否成立，不成立就一定不是素数，返回 false。

如果成立，且 $2\mid (p-1)$。

那么有 $(x^{\frac{p-1}{2}}{)}^2\equiv 1(\mathrm{mod}p)$。

这个式子长得和二次探测很像。

若 $x^{\frac{p-1}{2}}\mathrm{mod}p=1$，那么就可以继续判断，因为这样式子右边还是 $1$，还可以用二次探测。

若 $x^{\frac{p-1}{2}}\mathrm{mod}p=p-1$，那么无法判断 $x$ 是不是质数，返回 true。

其他情况，$x$ 一定不是质数，返回 false。

当 $2\nmid (p-1)$ 时，停止循环，返回 true。

虽然 Miller-Rabin 不一定正确，但可以多选取几个质数，这样正确率就很高了。

Miller-Rabin 代码如下：

ll qmul(ll a, ll b, ll p) { return (a * b - (ll)((long double)a / p * b) * p + p) % p; }

ll qpow(ll x, ll k, ll p) {
    x %= p;
    ll res = 1;
    while(k) {
        if(k & 1) res = qmul(res, x, p);
        x = qmul(x, x, p), k >>= 1;
    }
    return res;
}

const int pcnt = 12;
const ll p0[pcnt + 1] = { 0, 2, 3, 5, 7, 11, 13, 17, 19, 61, 2333, 4567, 24251 };

bool miller_rabin(ll x, ll p) {
    if(qpow(x, p - 1, p) != 1) return false;
    ll k = p - 1;
    while(!(k & 1)) {
        k >>= 1;
        ll t = qpow(x, k, p);
        if(t != 1 && t != p - 1) return false;
        if(t == p - 1) return true;
    }
    return true;
}

bool test(ll p) {
    if(p == 1) return false;
    lp(i, 1, pcnt) if(p == p0[i]) return true;
    lp(i, 1, pcnt) if(!miller_rabin(p0[i], p)) return false;
    return true;
}

Pollard-Rho

生日悖论

有 $n$ 个人，一年有 $d$ 天，存在两个人生日相同的概率是多少。

不难看出答案是 $p(n)=1-\frac{d!}{(d-n)!n^d}$，若 $d$ 为 365，那么 $n$ 为 $23$ 时 $p(n)$ 就超过了 $\frac{1}{2}$。

这启示我们如果选出若干个数，将它们两两作差，得到要分解的数 $n$ 的因数的概率很大。考虑 $n=p^2$ 时（最坏情况）。

设选出 $i$ 个数后使其存在两个数在模 $p$ 意义下相等的期望步数为 $E(i)$。

那么有
$$\begin{array}{rl}E(p)& =1\\ E(i)& =\frac{p-i}{p}(E(i+1)+1)+\frac{i}{p}\\ & =\frac{p-i}{p}E(i+1)+1\end{array}$$

$E(0)$ 大概是 $\sqrt{p}=n^{\frac{1}{4}}$，怎么得出来的呢？懒得推了，写个程序跑一遍就知道了。

但是这样做依旧不行，因为如果两两之间判断的话复杂度就又回到了 $O(n^{\frac{1}{2}})$，不是想要的结果。

伪随机数序列

Pollard-Rho 使用伪随机数生成序列（因为它有个特殊性质），假设现在生成了随机数 $x$，设 $f(x)=(x^2+C)\mathrm{mod}n$，那么 $x$, $f(x)$, $f(f(x))$… 构成了伪随机数序列。

它会形成一个形似 状的序列，即从某一点之后形成了循环。

图片来自 Pecco

它还有一个很好的性质。

若 $x\mid |a-b|$，那么 $f(a)-f(b)=a^2-b^2=(a+b)(a-b)$。

因此 $x\mid |f(a)-f(b)|$。

这个性质接下来会用到。

Pollard-Rho 算法

终于来了！

基于 Floyd 判环的 Pollard-Rho 算法

Pollard-Rho 通过不断将 $x$ 分解为两个因数的乘积，再将这两个因数继续分解，直到分解为质数。

每次执行 Pollard-Rho 时，先用 Miller-Rabin 判断要分解的数是否为质数，若是则说明找到了一个质因子。

Pollard-Rho 使用伪随机数序列，初始先随机出两个数 $t_1$ 和 $t_2$，令 $t_1\leftarrow f(t_2)$，即先让 $t_2$ 走一步。接下来每次循环 $t_1$ 走一步，$t_2$ 走两步。这样当 $t_1=t_2$ 时说明已经走完了一圈。

每次循环，判断是否 $gcd(|t_1-t_2|,x)\ge 2$，若是的话，说明找到了两个 $x$ 的因子，递归下去继续分解这两个因子。

由于之前说的性质，假设 $t_1$ 和 $t_2$ 之间距离为 $d$，判断了 $gcd(|t_1-t_2|,x)$，就相当于判断了序列上其他距离为 $d$ 的两个点，而 $t_1$ 和 $t_2$ 间的距离每次增加 $1$，因此相当于判断了序列上任意两个数。

另外若在一个环上没能分解 $x$，就再随机一次继续找下去。

这样做期望时间复杂度为 $O(n^{\frac{1}{4}}\log n)$。

void pollard_rho(ll x) {
    if(test(x)) {
        // 此时 x 为其中一个质因子，在这里做你想做的事
        return;   
    }
    while(true) {
        ll t1 = rd() % (x - 1) + 1, t2 = t1, b = rd() % (x - 1) + 1;
        t2 = (qmul(t2, t2, x) + t2 + b) % x;
        while(t1 != t2) {
            ll t = gcd(abs(t1 - t2), x);
            if(t != 1) {
                pollard_rho(t), pollard_rho(x / t);
                return;
            }
            t1 = (qmul(t1, t1, x) + t1 + b) % x;
            t2 = (qmul(t2, t2, x) + t2 + b) % x;
            t2 = (qmul(t2, t2, x) + t2 + b) % x;
        }
    }
}

倍增优化

我们可以将 $|t_1-t_2|$ 的结果乘起来，每隔 1，2，4，8…… 次循环求一次 gcd，因为如果 $\gcd (t,x)\ge 2$，那么有 $\gcd (kt\mathrm{mod}x)=\gcd ((k\mathrm{mod}x)(t\mathrm{mod}x),x)\ge 2$。

通常我们可以设一个最大距离，倍增超过这个距离后就一直保持这个距离。

特别的，当累乘的结果要变成 $0$ 时（因为是在模 $x$ 意义下的），就退出累乘，直接计算公因数。

期望时间复杂度 $O(n^{\frac{1}{4}})$ 。

void pollard_rho(ll x) {
    if(test(x)) return;
    while(true) {
        ll t1 = rd() % (x - 1) + 1, b = rd() % (x - 1) + 1, t2 = F(t1, b, x);
        int lim = 1;
        while(t1 != t2) {
            ll cnt = 1;
            lp(i, 1, lim) {
                ll tmp = qmul(cnt, abs(t1 - t2), x);
                if(!tmp) break;
                cnt = tmp;
                t1 = F(t1, b, x), t2 = F(F(t2, b, x), b, x);
            }
            ll t = gcd(cnt, x);
            if(t != 1) return pollard_rho(t), pollard_rho(x / t), void();
            lim = min(lim << 1, 128);
        }
    }
}

模板代码：

ll qmul(ll a, ll b, ll p) { return (a * b - (ll)((long double)a / p * b) * p + p) % p; }

ll qpow(ll x, ll k, ll p) {
    x %= p;
    ll res = 1;
    while(k) {
        if(k & 1) res = qmul(res, x, p);
        x = qmul(x, x, p), k >>= 1;
    }
    return res;
}

const int pcnt = 12;
const ll p0[pcnt + 1] = { 0, 2, 3, 5, 7, 11, 13, 17, 19, 61, 2333, 4567, 24251 };

bool miller_rabin(ll x, ll p) {
    if(qpow(x, p - 1, p) != 1) return false;
    ll k = p - 1;
    while(!(k & 1)) {
        k >>= 1;
        ll t = qpow(x, k, p);
        if(t != 1 && t != p - 1) return false;
        if(t == p - 1) return true;
    }
    return true;
}

bool test(ll p) {
    if(p == 1) return false;
    lp(i, 1, pcnt) if(p == p0[i]) return true;
    lp(i, 1, pcnt) if(!miller_rabin(p0[i], p)) return false;
    return true;
}

ll gcd(ll a, ll b) { return a ? gcd(b % a, a) : b; }

mt19937 rd(time(NULL));

ll F(ll x, ll b, ll p) { return (qmul(x, x, p) + x + b) % p; }

void pollard_rho(ll x) {
    if(test(x)) return;
    while(true) {
        ll t1 = rd() % (x - 1) + 1, b = rd() % (x - 1) + 1, t2 = F(t1, b, x);
        int lim = 1;
        while(t1 != t2) {
            ll cnt = 1;
            lp(i, 1, lim) {
                ll tmp = qmul(cnt, abs(t1 - t2), x);
                if(!tmp) break;
                cnt = tmp;
                t1 = F(t1, b, x), t2 = F(F(t2, b, x), b, x);
            }
            ll t = gcd(cnt, x);
            if(t != 1) return pollard_rho(t), pollard_rho(x / t), void();
            lim = min(lim << 1, 128);
        }
    }
}

模板题：洛谷 [P4718]

​