Miller Rabin素数检测与Pollard Rho算法
一些前置知识可以看一下我的联赛前数学知识
如何判断一个数是否为质数
方法一：试除法
扫描\(2\sim \sqrt{n}\)之间的所有整数，依次检查它们能否整除\(n\)，若都不能整除，则\(n\)是质数，否则\(n\)是合数。
代码
bool is_prime(int n){
if(n<2) return 0;
int m=sqrt(n);
for(int i=2;i<=m;i++){
if(n%i==0) return 0;
}
return 1;
}
方法二、线性筛
用 \(O(n)\) 的复杂度筛出所有的质数，然后 \(O(1)\) 判断
代码
const int maxn=1e6+5;
int pri[maxn];
bool not_prime[maxn];
void xxs(int n){
not_prime[0]=not_prime[1]=1;
for(int i=2;i<=n;++i){
if(!not_prime[i]){
pri[++pri[0]]=i;
}
for (int j=1;j<=pri[0] && i*pri[j]<=n;j++){
not_prime[i*pri[j]]=1;
if(i%pri[j]== 0) break;
}
}
}
但是，我们来看一下这道题 LOJ#143. 质数判定
在 \(5s\) 的时间内判断 \(10^5\) 个 \(10^{18}\) 级别的数是否是质数
以上两种方法显然都不可行
所以我们要用到一种更高效的算法 \(Miller\ Rabin\)
Miller Rabin素数检测
根据费马小定理
若 \(p\) 为素数，对于任意整数 \(a\)，都有
\(a^{p-1} \equiv 1(mod\ p)\)
它的逆定理在大多数情况下是成立的
所以，可以每一次随机挑选一些数 \(a\)
判断是否存在 \(a^{p-1} \equiv 1(mod\ p)\)
只要有一个 \(a\) 不满足条件，就将其标记为合数
否则标记为质数
代码
#include<cstdio>
#include<ctime>
#include<cstdlib>
#define rg register
typedef long long ll;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs>>=1LL;
}
return nans;
}
int main(){
rg ll aa,bb;
rg bool jud=1;
while(scanf("%lld",&aa)!=EOF){
jud=1;
if(aa==1){
printf("N\n");
continue;
}
for (rg int i=1;i<=100;i++) {
bb=1LL*rand()*rand()%(aa-1)+1;
if(ksm(bb,aa-1,aa)!=1){
jud=0;
break;
}
}
if(jud) printf("Y\n");
else printf("N\n");
}
return 0;
}
但是也有例外，即存在一种极端反例卡迈克尔数（一种合数）
对于任何卡迈克尔数，费马定理都成立
虽然这种数很少，但是还是有可能会被卡（只有\(33\)分）
所以我们需要借助另一个定理来优化它
二次探测定理：若 \(p\) 为质数且 \(x∈(0,p)\)，则方程 \(x^2≡1(mod\ p)\) 的解为 \(x_1=1,x_2=p-1\)
对于任意一个奇素数 \(p\) ，\(p-1\) 一定可以写成 \(r2^t\) 的形式
因此我们可以把费马小定理写成 \(a^{r2^t}\equiv 1(mod\ p)\) 的形式
一开始，我们先求出 \(a^r mod\ p\) 的值
然后每一次给这个值平方，一共平方 \(t\) 次
算一下每次得出来的结果是否满足二次探测定理
如果不满足，说明这个数不是质数
最后再看一下最终的值是否满足费马小定理即可
时间复杂度 \(klogn\)
其中 \(k\) 为每次检测的次数
事实证明，这个算法出错的概率非常小
为 \(\frac{1}{4^k}\)
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<iostream>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
typedef long long ll;
const int times=7;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs>>=1LL;
}
return nans;
}
bool check(rg ll a,rg ll r,rg ll t,rg ll mod){
ll nans=ksm(a,r,mod),tmp=nans;
for(rg ll i=1;i<=t;i++){
nans=gsc(tmp,tmp,mod);
if(nans==1 && tmp!=1 && tmp!=mod-1) return 0;
tmp=nans;
}
if(tmp==1) return 1;
else return 0;
}
bool Miller_Robin(rg ll n){
if(n==2) return 1;
if(n<2 || (n&1LL)==0) return 0;
rg ll t=0,r=n-1;
while((r&1LL)==0){
r>>=1;
t++;
}
for(rg int i=1;i<=times;i++){
rg ll a=rand()%(n-1)+1;
if(!check(a,r,t,n)) return 0;
}
return 1;
}
ll x;
int main(){
srand(time(0));
while(scanf("%lld",&x)!=EOF){
if(Miller_Robin(x)) printf("Y\n");
else printf("N\n");
}
return 0;
}
Pollard Rho算法
这道题不仅要判断是否是质数，还要求输出最大质因子
判质数的话用 \(Miller\ Rabin\) 判一下就好了
关键是怎么找出最大质因子
其实还是随机的思想
我们每一次随机一个数 \(n\)，判断 \(n\) 是不是 \(p\) 的因数
如果是，就分成两个子问题 \(n\) 和 \(\frac{p}{n}\) 递归求解
如果当前的数为质数就更新答案
但是这样随机概率非常小
利用生日悖论，采用组合随机采样的方法,满足答案的组合比单个个体要多一些.这样可以提高正确率
具体方法是随机两个整数 \(n\)，\(m\)
判断 \(|n-m|\) 是不是 \(p\) 的因数
看起来没有什么不同，实际上效率提高了不少
剩下的就在于怎么构造随机的 \(n,m\) 了
构造的方法会影响到我们程序的效率
实践证明，直接 \(rand\) 效率极低
而构造一种形如 \(f(x)=x^2+c\) 的数列效率很高
我们首先 \(rand\) 一个 \(x_1\)
然后令 \(f(x_2)=x_1^2+c\)
\(f(x_3)=x_2^2+c\)
依次类推，每次把序列相邻两项作差判断是否是 \(p\) 的因数
但是因为我们是在模意义下进行运算，所以这个序列一定有循环节
所以当遇到循环节的时候我们就退出循环，重新构造一个数列
因为每次求 \(gcd\) 的开销太大
所以我们可以先把相邻两项的差连乘起来，这是不影响结果的
当累乘的次数等于我们设定的一个值时再进行求 \(gcd\) 的运算
一般把这个值设为 \(2^k\)
代码
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<cmath>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
typedef long long ll;
const int times=10;
int nt;
ll x,ans;
ll gsc(rg ll ds,rg ll zs,rg ll mod){
return ((unsigned long long)(ds*zs)-(unsigned long long)((long double)ds/mod*zs)*mod+mod)%mod;
}
ll ksm(rg ll ds,rg ll zs,rg ll mod){
rg ll nans=1;
while(zs){
if(zs&1LL) nans=gsc(nans,ds,mod);
ds=gsc(ds,ds,mod);
zs>>=1LL;
}
return nans;
}
bool check(rg ll a,rg ll r,rg ll t,rg ll mod){
ll nans=ksm(a,r,mod),tmp=nans;
for(rg ll i=1;i<=t;i++){
nans=gsc(tmp,tmp,mod);
if(nans==1 && tmp!=1 && tmp!=mod-1) return 0;
tmp=nans;
}
if(tmp==1) return 1;
else return 0;
}
bool Miller_Robin(rg ll n){
if(n==2) return 1;
if(n<2 || (n&1LL)==0) return 0;
rg ll t=0,r=n-1;
while((r&1LL)==0){
r>>=1;
t++;
}
for(rg int i=1;i<=times;i++){
rg ll a=rand()%(n-1)+1;
if(!check(a,r,t,n)) return 0;
}
return 1;
}
ll gcd(rg ll aa,rg ll bb){
if(bb==0) return aa;
return gcd(bb,aa%bb);
}
ll rp(rg ll x,rg ll y){
rg int t=0,k=1;
rg ll v0=rand()%(x-1)+1,v=v0,d,s=1;
while(1){
v=(gsc(v,v,x)+y)%x;
s=gsc(s,std::abs(v-v0),x);
if(v==v0 || !s) return x;
if(++t==k){
d=gcd(s,x);
if(d!=1) return d;
s=1,v0=v,k<<=1;
}
}
}
void dfs(rg ll n){
if(n<=ans || n==1) return;
if(Miller_Robin(n)){
ans=n;
return;
}
rg ll now=n;
while(now==n) now=rp(n,rand()%n+1);
while(n%now==0) n/=now;
dfs(n),dfs(now);
}
ll solve(rg ll n){
ans=1;
dfs(n);
return ans;
}
int main(){
srand(time(0));
scanf("%d",&nt);
while(nt--){
scanf("%lld",&x);
if(Miller_Robin(x)) printf("Prime\n");
else printf("%lld\n",solve(x));
}
return 0;
}
