zoukankan      html  css  js  c++  java
  • Miller_Rabin 素数判定算法

    我们先看一道模板题
    LOJ #143 素数判定
    也许你会想,这不就是素数判定吗,直接上试除法,就这!
    然而,这道题目有\(10^5\)次询问,每次询问的\(x \le 10^{18}\)
    普通的\(O(\sqrt x)\)的试除法根本过不去,这是就轮到我们的主角:Miller_Rabin算法登场了

    2个重要定理

    1.费马小定理
    对于任意的质数\(p\),满足\(a^{p-1} \equiv 1(mod\ p)\)

    2.二次探测定理
    对于任意的质数\(p\)及整数\(x\),若满足\(x^2 \equiv 1(mod\ p)\),则一定有\(x \equiv 1或p-1(mod\ p)\)

    算法流程

    想必你一定会想用以上2个定理的逆命题来判定质数吧。
    遗憾的是,二次探测定理逆命题的反例非常多,而对于费马小定理,也已被证明存在无数多个合数满足费马小定理逆命题,叫做\(Carmichael\)数。
    但好消息是,\(Carmichael\)数的密度极小,\(1e8\)范围内只有\(255\)个,并且在进行次数较多的二次探测定理逆否命题判定后,\(Carmichael\)数都会被判定出来。

    于是我们考虑通过结合二者来进行判定:
    令我们需要判定的\(x=2^s*t+1\)
    然后通过多次寻找一个较小的数\(a\),先计算出\(a^t\ mod\ x\),在进行\(s\)次平方,每次平方的同时使用二次探测定理进行判定
    平方\(s\)次后,最终得到的数\(=(a^t)^{2^s}=2^{x-1}\),于是我们用费马小定理的逆否命题来进行判定

    容易看出,当\(a\)是较小的质数时,判定出错的概率会更小
    因为有时间复杂度的限制,所以我们选择使用最小的10个质数作为\(a\)依次进行上述判定,顺带的,我们可以先将这些质数的倍数特判以降低一些常数
    这样一来,我们进行了10次费马小定理的判定与数十次二次探测定理,这些结合起来之后,算法的错误概率几乎已经可以忽略不计了,因此\(Miller\_Rubin\)是一个可靠的算法

    需要注意的是,因为\(x\)\(long\ long\)范围的数,因此在进行快速幂时,直接进行乘法可能会爆\(long\ long\)
    所以我们选择将乘法用类似快速幂的方式,使用\(log\)次加法来实现
    于是我们就得到了一个十分优秀,足以通过本题的素数判定算法

    code

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    ll x;
    ll pr[20]={2,3,5,7,11,13,17,19,23,29};
    inline ll add(ll x,ll y,ll mod){return (x+y>=mod)?x+y-mod:x+y;}
    inline ll mul(ll a,ll b,ll mod){
    	ll ret=0;
    	for(;b;a=add(a,a,mod),b>>=1) if(b&1) ret=add(ret,a,mod);
    	return ret;
    }
    inline ll fsp(ll a,ll b,ll mod){
    	ll ret=1;
    	for(;b;a=mul(a,a,mod),b>>=1) if(b&1) ret=mul(ret,a,mod);
    	return ret;
    }
    inline bool Miller_Rabin(ll x){
    	if(x<2) return false;
    	for(int i=0;i<10;++i){
    		if(x==pr[i]) return true;
    		if(x%pr[i]==0) return false;
    	}
    	int s=0;ll t=x-1;
    	while(!(t&1)){s++;t>>=1;}
    	for(int i=0;i<10;++i){
    		ll a=pr[i];
    		ll cur=fsp(a,t,x),nxt;
    		for(int j=0;j<s;++j,cur=nxt){
    			nxt=mul(cur,cur,x);
    			if(cur!=1&&cur!=x-1&&nxt==1) return false;
    		}
    		if(cur!=1) return false;
    	}
    	return true;
    }
    int main(){
    	while(scanf("%lld",&x)!=EOF){
    		if(Miller_Rabin(x)) puts("Y");
    		else puts("N");
    	}
    	return 0;
    }
    
  • 相关阅读:
    关于求 p_i != i and p_i != i+1 的方案数的思考过程
    poj 3041 Asteroids 二分图最小覆盖点
    poj 1325 Machine Schedule 最小顶点覆盖
    poj 1011 Sticks 减枝搜索
    poj 1469 COURSES 最大匹配
    zoj 1516 Uncle Tom's Inherited Land 最大独立边集合(最大匹配)
    Path Cover (路径覆盖)
    hdu 3530 SubSequence TwoPoint单调队列维护最值
    zoj 1654 Place the Rebots 最大独立集转换成二分图最大独立边(最大匹配)
    poj 1466 Girls and Boys 二分图最大独立子集
  • 原文地址:https://www.cnblogs.com/tqxboomzero/p/14038921.html
Copyright © 2011-2022 走看看