【旧文/OI】素数判定与质因子寻找（PR）算法

发表于： 更新于： 分类于：
字数： 1125 阅读：≈ 3分钟 浏览：

本博客源自我的洛谷博客 https://www.luogu.com.cn/article/jbqa6q5t 和 https://www.luogu.com.cn/article/bo0gb1ur （现在改成专栏了），写作时间 2021-07-30 21:59/2021-08-13 20:08

OI-wiki上关于素数的介绍

matrix67大佬对素数的介绍

素数是个啥

不会有人不知道吧qaq

素数 就是大于一的、因数只有一和它自身的数，例如：

2357111317192329……

判定素数

我不认为我能讲的比 OI-wiki 上的更易懂，在这里精炼的记一下公式：

枚举法

$\Theta(\sqrt n)$ 复杂度，这个应该都熟知了吧：

1
2
3
4
5

bool isp(int p){
    for(int i=2;i*i<=p;i++)
        if(p%i==0)	return 0;
    return 1;
}

Miller-Rabin 素数测试

这个是 OI 中很常用的素性测试

费马小定理：

$p$ 是素数 $\implies a^p \equiv a\pmod p$

二次探测定理：

$p$是素数 & $ x^2\equiv1\pmod p\implies (\ x\equiv1\pmod p\ \operatorname{or}\ x\equiv p-1\pmod p\ )$

Miller Rabin

若 $p$ 是素数，$\gcd(a,p)=1$：

根据费马小定理：$a^{p-1}\equiv1\pmod p$

使用二次探测定理：$a^\frac{p-1}{2}\equiv1\ \operatorname{or}\ p-1\pmod p$

若 $\frac{p-1}{2}$ 是偶数，且 $a^\frac{p-1}{2}\equiv1\pmod p$，则可继续使用二次探测定理。

根据这个原理，就可以写出 Miller-Rabin 的板子啦！

奉上我丑陋的代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

bool isp(ll p){//based on Miller-Rabin
    if(p<3||p%2==0)	return p==2;
    for(int i=1;i<=20;i++){ //repeat 20 times
        ll a=rand()%(p-2)+2,x=p-1;bool flag=1,exf=1;
        if(p%a==0)	return 0;
        if(pow(a,x,p)!=1ll)return 0;x/=2;
        for(ll r;exf;x/=2ll){
            r=pow(a,x,p);if(x%2ll)exf=0;
            if(r!=1){flag=(r==p-1);break;}
        }if(!flag)return 0;
    }return 1;
}

pow(a,b,mod) 即带 取模的快速幂 ，注意其中乘法应使用带取模的快速乘。

Lucas–Lehmer 算法（针对梅森素数的素数判定）

梅森素数是形如 $2^k-1$ 的素数，这个算法能判断 $2^k-1$ 是不是素数。

首先构造一个序列： $$L_i=\begin{cases}4,&i=0 \\ L_{i-1}^2-2,&\text{otherwise.}\end{cases}$$

则 $2^k-1$ 是素数当且仅当 $L_{p-2}\equiv0\pmod{2^k-1}$

故可写出代码：

1
2
3
4
5
6

bool isp(int k){ //test on (2^k)-1
    ll p=(1ll<<k)-1ll,Li=4ll;
    for(int i=1;i<=k-2;i++)
        Li=(mul(Li,Li,p)-2+p)%p;
    return Li==0;
}

其中 mul(a,b,p) 是快速乘，见 这里 。

参考： https://baike.baidu.com/item/卢卡斯-莱默检验法/

https://www.wx1986.com/program/867.html/

AKS 算法

AKS 算法也是一种 “一般的、多项式的、确定性的和无仰赖的” 判定素数的方法，复杂度比 MR 差，但胜在不是随机的（其实我觉得也没必要，MR 多循环几次一般就够了），我不懂原理，所以链接：

https://zhuanlan.zhihu.com/p/346563055/

Pollard-Rho算法

参考： https://zhuanlan.zhihu.com/p/267884783/

作者讲的非常清楚！

给个Floyd判环法的demo： https://visualgo.net/zh/cyclefinding/

OI-wiki的介绍： https://oi-wiki.org/math/number-theory/pollard-rho/

主体思路

Pollard-Rho 算法是一个随机算法，主体思路大概是……：

假设我要找 $a$ 的质因数，随机一个数 $p$ 判断它与 $a$ 是否互质，如果不是就递归找 $\gcd(a,p)$ 和 $\dfrac a {\gcd(a,p)}$ 的因数，否则重找。

然后 Pollard-Rho 使用独特的随机函数，提升更有可能不与 $a$ 互质的 $p$ 的出现概率，且尽量让 $p$ 不重复。

就是一个看脸的算法啦~，详细实现建议去上面“ 参考 ”里看看。

想通了不是很难，代码也不是很长啦。

代码实现

附上我丑陋的代码：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

ll PR(ll N){	ll c,t,r,pd,tmp;      //refer to: https://zhuanlan.zhihu.com/p/267884783
    #define rho(a) (mul(a,a,N)+c)%N   //a special random number generator by Pollard,
    if(isp(N))return N;if(N==4)return 2;//if N==4 no matter what c, no solution
    while(1){                         //try this many times
        c=rand()%(N-2)+1;        //Floyd's way to determine loop (simplified). demo on https://visualgo.net/zh/cyclefinding
        t=rho(0);r=rho(rho(0));pd=abs(t-r);      //t for turtle, r for rabbit, pd to gather products, rho(pre) to get the next number
        for(int lim=1;t!=r;lim=min(128,lim<<1)){   //using this to get every number with distant of d on the loop
            for(int i=0;i<lim;i++){                //set const distant to 128
                t=rho(t);r=rho(rho(r));tmp=mul(pd,abs(t-r),N);
                if(t==r||tmp==0)break;pd=tmp;      //use pd to gather products to reduce #gcd used times
            }
            tmp=gcd(pd,N);if(tmp>1){return max(PR(tmp),PR(N/tmp));}
        }
    }
}

isp(N) 是用 Miller Rabin 判断指数的板子，详见：

mul(a,b,p) 是快速乘，见 https://oi-wiki.org/math/number-theory/mod-arithmetic/#快速乘

大概……就这样