zoukankan      html  css  js  c++  java
  • 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;
    }
    
  • 相关阅读:
    mysql中timestamp自动更新
    Jquery JQZoom Plugin 放大鏡效果 Two
    ASP.NET MVC 之路(一)HtmlHelper 擴展控件
    我的第一个控件
    Visifire Silverlight Chart 比 Silverlight Toolkit Chart 好的兩點
    Jquery Magnify Plugin 放大鏡效果 One
    Asp.Net MVC 之路(二)Html Helper擴展之二
    一个信仰从此在这里开始(仅此篇是非技术)
    HelloWord之WCF文件的猫咪形式
    第一话 Asp.Net MVC 3.0【Hello World!】
  • 原文地址:https://www.cnblogs.com/liuchanglc/p/14245195.html
Copyright © 2011-2022 走看看