LGP4718 [LG TPLT] Pollard-Rho 学习笔记

LGP4718 [LG TPLT] Pollard-Rho 学习笔记

Luogu Link

题意简述

\(T\) 组询问，每次询问给定一个整数 \(n\)。若 \(n\) 为质数则报告之，否则求其最大的质因子。

\(T\le 350,2\le n\le 10^{18}\)。

做法解析

更快的素性测试？

啧，要是按 \(\sqrt{n}\) 的试除法，我们似乎连 \(n\) 是不是质数都判不掉啊，更不谈质因数分解了。所以我们需要使用Miller-Rabin素性测试来解决此题。

什么是Miller-Rabin素性测试？

费马素性测试

首先我们要知道费马小定理：若 \(p\) 为素数，则 \(\forall 1\le a\le p\)，有 \(a^{p-1}\equiv 1\pmod{p}\)。

那么它有逆否命题：若 \(\exists 1\le a\le p\)，有 \(a^{p-1}\not\equiv 1\pmod{p}\)。

那么我们可以在 \([1,p-1)\) 里随便选一些数并快速幂之。如果存在 \(a^{p-1}\not\equiv 1\pmod{p}\) 则 \(p\) 一定为合数，否则 \(p\) 大概率为素数。这被称为费马素性检验。

不幸地，存在一些合数 \(p\) 也满足 \(a^{p-1}\not\equiv 1\pmod p\)。而且这类数的数量不支持打表判掉它们。所以我们得再验证一道。

二次探测定理与Miller-Rabin

我们有二次探测定理：若 \(p\) 为奇素数，则 \(x^2\equiv 1\pmod{p}\) 的解为 \(x\equiv 1\pmod{p}\) 或 \(x\equiv p-1\pmod{p}\)。

证明：移项上柿可得 \(x^2-1\equiv 0\pmod{p}\)。再平方差公式一下即可。

所以我们可以这么优化算法（我们默认我们要检验的 \(p\) 是奇数，因为偶数可以直接判掉）：

• 我们选取若干个底数 \(a\)，对于当前这个 \(a\)，我们设 \(d=p-1\)，则我们快速幂判断是否有 \(a^d\equiv 1\pmod{p}\)，如果不是，则 \(p\) 为合数。
• 因为 \(p\) 为奇数，所以 \(d=p-1\) 必为偶数，我们可以把刚才的柿子看作 \((a^{\frac{d}{2}})^2\equiv 1\pmod{p}\)。
• 那么我们就来检验一下 \(a^{\frac{d}{2}}\) 是否是 \(\pm 1\)，如果否，则 \(p\) 必为合数，如果是 \(1\) 且 \(d\) 不为奇数，那这个判断过程还能递归下去。

在随机选取 \(k\) 个底数的情况下错误率期望是 \(\frac{1}{4^k}\) 的。不过这是OI，数的值域不太会干过 \(2^{64}\)。按照上述方法，选择前十二个质数作为底数，在 \(2^{78}\) 内是确定性的。

另外，如果值域只有 int 级别，选择前四个质数即可，如果到了 uint，则需要选择前五个质数。

Pollard-Rho

好了终于到正题了。\(O(\sqrt{n})\) 试除法太慢了，我们得优化。

Pollard-Rho是一种随机算法，其目标就是快速用随机数试出 \(n\) 的因子。接下来我们默认 \(n\) 是合数（我们在跑Pollard-Rho之前先用Miller-Rabin把 \(n\) 是素数的情况直接判掉）。

我们知道一个合数的最小质因子 \(m\) 满足 \(m\le \sqrt{n}\)。生日悖论又告诉我们，对于一个理想的取值为 \([1,n]\) 的随机数生成器，生成 \(\sqrt{\frac{\pi n}{2}}\) 个数可以期望得到两个数相同。

所以我们如果随机生成 \(\sqrt[4]{n}\) 个数，期望就存在 \(a,b\) 满足 \(a\equiv b\pmod{m}\)。那么 \(|a-b|\) 就是 \(p\) 的倍数，我们就可以通过求 \(\gcd(|a-b|,n)\) 得到一个因子。

问题是我们要枚举所有对 \((a,b)\)，复杂度又回到 \(O(\sqrt{n})\) 了，又因为 \(\gcd\)，在此基础上还多了个 \(\log n\)。考虑优化。

优秀的函数和正确的期望

Pollard-Rho 设计了这样一种函数：\(a_i=f(a_{i-1})=a_{i-1}^2+c\)。这个函数有一个特点：若 \(a_i=a_j\)，则 \(a_{i+1}=f(a_i)=f(a_j)=a_{j+1}\)。实际上，这个数列形成一种混循环的状态，出现上述情况就是进环了（出现混循环也是非常自然的，因为在模意义下值域是有限的，所以总有走到之前走过的值域的那一刻）。

Pollard-Rho接下来是这么实现的：第 \(k\) 次操作让一个数为 \(a_k\)，另一个 \(a_{2k}\)，如果 \(a_k=a_{2k}\) 那就毫无疑问进入环了，操作跳出（显然这个过程肯定是有尽头的，因为 \(k\) 慢慢增长一定会增长到环长的若干倍），否则就检查 \(\gcd(|a_k-a_{2k}|)\)。

如果跑了一次 \(\rho\) 形的判断还没找出不等于 \(1\) 的 \(\gcd\)，就重新随机一个 \(a_0\) 和 \(c\)。

期望上，上述混循环的链长和环长都是 \(\sqrt{m}\) 的（根据生日悖论），检查的数对也可以认为是 \(O(\sqrt{m})\) 级别（所以说，“若 \(a_i=a_j\)，则 \(a_{i+1}=a_{j+1}\)”也可以理解为：一种环上的距离对应一个结果。有 \(\sqrt{m}\) 种有意义的距离，加上环外我们就检查了 \(O(\sqrt{m})\) 次 \(\gcd\)，期望上就是对的）。

Action开卡

然而 \(\gcd\) 做的太多，常熟就太大了。

在这种 \(\gcd(x_i,n)\) 大多数都是 \(1\)，只有少数 \(\gcd(x_i,n)\) 不为 \(1\)，我们只关心 \(\gcd(x_i,n)\) 的非一值的积的情况下，我们不妨做 \(\gcd(\prod x_i,n)\)，这么做的唯一问题在于如果 \(\prod x_i | n\) 了就糟糕了，我们能分解出来的 \(\gcd\) 就被吞了。所以我们的策略是以 \(\log n\) 作为 \(\prod x_i\) 的段长，这样既平衡了 \(\gcd\) 的复杂度，又让 \(\gcd\) 被吞的可能性不至于太高。

（另外，实现时步长是从 \(1\) 慢慢加的，似乎可以提高效率？）

代码实现

终于到代码了。
polaro 函数中把 \(n=4\) 判掉是因为我们随机 \(a_k\) 和 \(c\) 的值域范围是 \([3,n-1]\)，\(n=4\) 的时候这个值域就是 \([3,3]\)，没有随机性了。

注意这个 fgcd，它通过快速约掉那些潜在含有的 \(2\)，比普通的更相减损要快不少（辗转相除虽然复杂度 \(O(\log n)\)，但是对大数字取模太慢），也有 \(O(\log n)\) 的复杂度保底（更相减损最坏 \(O(n)\)）。虽然说，CPP自己实现的求 \(\gcd\) 就是这样的，但是现在赛场上只有CPP14，所以还是得准备自己写。

#include <bits/stdc++.h>
using namespace std;
using namespace obasic;
namespace omathe{
    const int mrps[12]={2,3,5,7,11,13,17,19,23,29,31,37};
	lolo fgcd(lolo a,lolo b){
        if(!a||!b)return a|b;
        int t=binctzl(a|b);b>>=t;
        for(;b;a-=b){
            a>>=binctzl(a);
            if(a<b)swap(a,b);
        }
        return a<<t;
    }
    molo fpow(molo a,molo b,molo p){
        a%=p;molo res=1;
        for(;b;(a*=a)%=p,b>>=1)if(b&1)(res*=a)%=p;
        return res;
    }
    bool mrchk(lolo p,int a){
        lolo d=p-1,tmp=fpow(a,d,p);
        if(tmp!=1)return false;
        for(;!(d&1);){
            d>>=1,tmp=fpow(a,d,p);
            if(tmp==p-1)return true;
            if(tmp!=1)return false;
        }
        return true;
    }
    bool milrab(lolo p){
        if(p<3||!(p&1))return p==2;
        if(p>40){for(auto a : mrps)if(!mrchk(p,a))return 0;return 1;}
        for(auto a : mrps)if(a==p)return 1;return 0;
    };
    molo Plrf(molo x,molo c,molo p){return (x*x+c)%p;}
    lolo polaro(lolo n){
        if(n==4)return 2;
        lolo d=0;lolouid cuid(3,n-1);
        lolo l=cuid(mrnd),r=l,c=cuid(mrnd);
        l=Plrf(l,c,n);r=Plrf(l,c,n);
        for(int lim=1;l!=r;lim<<=1,minner(lim,(int)log2(n)+1)){
            molo cpd=1;
            for(int i=0;i<lim;i++){
                (cpd*=abs(l-r))%=n;if(!cpd)break;
                l=Plrf(l,c,n),r=Plrf(Plrf(r,c,n),c,n);
            }
            d=fgcd(cpd,n);if(d!=1)return d;
        };
        return n;
    }
    lolo polaros(lolo n){lolo res=n;for(;res==n;res=polaro(n));return res;}
};
using namespace omathe;
lolo solve(lolo n){
    if(milrab(n))return n;
    lolo d1=polaros(n),d2=n/d1,res=0;
    maxxer(res,max(solve(d1),solve(d2)));
    return res;
}
lolo N,ans;
void mian(){
    readi(N),ans=solve(N);
    if(N==ans)puts("Prime");
    else writil(ans);
}
int Tcn;
int main(){
    readi(Tcn);
    while(Tcn--)mian();
    return 0;
}

后记

我是大常熟选手.jpg。

后记的后记

两个月后，我终于找到了常熟大的罪魁祸首：

return (x*x%p+c)%p;
//点击输入文本
//双击修改副标题