本博客源自我的洛谷博客https://www.luogu.com.cn/article/jbqa6q5t和https://www.luogu.com.cn/article/bo0gb1ur (现在改成专栏了),写作时间 2021-07-30 21:59/2021-08-13 20:08
素数是个啥
不会有人不知道吧qaq
素数就是大于一的、因数只有一和它自身的数,例如:
| 2 | 3 | 5 | 7 | 11 | 13 | 17 | 19 | 23 | 29 | …… |
|---|
判定素数
我不认为我能讲的比 OI-wiki 上的更易懂,在这里精炼的记一下公式:
枚举法
\Theta(\sqrt n) 复杂度,这个应该都熟知了吧:
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 的板子啦!
奉上我丑陋的代码:
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}
故可写出代码:
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://baike.baidu.com/item/AKS素数测试/22735657/
https://zhuanlan.zhihu.com/p/346563055/
参考: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 不重复。
就是一个看脸的算法啦~,详细实现建议去上面“参考”里看看。
想通了不是很难,代码也不是很长啦。
附上我丑陋的代码:
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/quick-pow/#_10/
大概……就这样,如果有问题,可以到评论中问我噢!(/虽然我可能看不到/)