Miller Rabin素数检测与Pollard Rho算法

一些前置知识可以看一下我的联赛前数学知识

如何判断一个数是否为质数

方法一:试除法

扫描(2sim 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
");
			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
");
		else printf("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
");
		else printf("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
");
		else printf("%lld
",solve(x));
	}
	return 0;
}
原文地址:https://www.cnblogs.com/liuchanglc/p/14245195.html